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A New Method for Recognizing Quadric Surfaces from Range Data and 
Its Application to Telerobotics and Automation 

(phase II) 

by 

Nicolas Alvertos* and Ivan D’Cunha** 

Abstract 

The problem of recognizing and positioning of objects in three-dimensional space 
is important for robotics and navigation applications. In recent years, digital range 
data, also referred to as range images or depth maps, have been available for the 
analysis of three-dimensional objects owing to the development of several active range 
finding techniques. The distinct advantage of range images is the explicitness of the 
surface information available. Many industrial and navigational robotics tasks will be 
more easily accomplished if such explicit information can be efficiently interpreted. 

In this research, a new technique based on analytic geometry for the recognition 
and description of three-dimensional quadric surfaces from range images is pressented. 
Beginning with the explicit representation of quadrics, a set of ten coefficients are 
determined for various three-dimensional surfaces. For each quadric surface, a unique 
set of two-dimensional curves which serve as a feature set is obtained from the various 
angles at which the object is intersected with a plane. Based on a discriminant 
method, each of the curves is classified as a parabola, circle, ellipse, hyperbola, or a 
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line. Each quadric surface is shown to be uniquely characterized by a set of these 
two-dimensional curves, thus allowing discrimination from the others. 

Before the recognition process can be implemented, the range data have to 
undergo a set of pre-processing operations, thereby making it more presentable to 
classification algorithms. One such pre-processing step is to study the effect of median 
filtering on raw range images. Utilizing a variety of surface curvature techniques, reli- 
able sets of image data that approximate the shape of a quadric surface are determined. 
Since the initial orientation of the surfaces is unknown, a new technique is developed 
wherein all the rotation parameters are determined and subsequently eliminated. This 
approach enables us to position the quadric surfaces in a desired coordinate system. 

Experiments were conducted on raw range images of spheres, cylinders, and 
cones. Experiments were also performed on simulated data for surfaces such as hyper- 
boloids of one and two sheets, elliptical and hyperbolic paraboloids, elliptical and 
hyperbolic cylinders, ellipsoids and the quadric cones. Both the real and simulated 
data yielded excellent results. Our approach is found to be more accurate and compu- 
tationally inexpensive as compared to traditional approaches, such as the three- 
dimensional discriminant approach which involves evaluation of the rank of a matrix. 

Finally, we have proposed one other new approach, which involves the formula- 
tion of a mapping between the explicit and implicit forms of representing quadric sur- 
faces. This approach, when fully realized, will yield a three-dimensional discriminant, 
which will recognize quadric surfaces based upon their component surface patches. 
This approach is faster than prior approaches and at the same time is invariant to pose 
and orientation of the surfaces in three-dimensional space. 


Table of Contents 


LIST OF TABLES 

LIST OF FIGURES 

CHAPTER 1: INTRODUCTION 

1.1 Introduction 

1.2 Range Image and Data Acquisition 

1.3 Definition of Object Recognition Problem 

1.4 Objectives and Organization of the Dissertation 

CHAPTER 2: BACKGROUND 

2.1 Introduction 

2.2 Literature Review 

2.3 Differential Geometry of Surfaces: Mean and Gaussian 

Curvatures 

2.4 Three-Dimensional Discriminant 

2.4. 1 Translation 

2.4.2 Rotation 

CHAPTER 3: QUADRIC SURFACE REPRESENTATION 

3.1 Introduction 

3.2 Quadric Surface Description 


vi 

ix 

3 

1 

1 

2 

3 

5 

5 

5 

8 

12 

13 

13 

19 

19 

19 


m 


3.3 Recognition Scheme 


23 


3.3.1 Median Filtering 

3.3.2 Segmentation 

3.3.3 Three-Dimensional Coefficients Evaluation 

3.3.4 Evaluation of the Rotation Matrix 

3.3.5 Product Terms Elimination Method 

3.3.6 Translation of the Rotated Object 

CHAPTER 4: QUADRIC SURFACE CHARACTERIZATION AND 
RECOGNITION 

4.1 Introduction 

4.2 Two-Dimensional Discriminant 

4.3 Quadric Surface Description and Representation 

4.3.1 Ellipsoid 

4.3.2 Circular (Elliptic) cylinder 

4.3.3 Sphere 

4.3.4 Quadric circular (elliptic) cone 

4.3.5 Hyperboloid of one sheet 

4.3.6 Hyperboloid of two sheets 

4.3.7 Elliptic paraboloid 

4.3.8 Hyperbolic paraboloid 

4.3.9 Hyperbolic cylinder 

4.3.10 Parabolic cylinder 

4.3.11 Parallelepiped 


24 

28 

28 

33 

38 

44 


47 

47 

48 

49 
49 
55 
62 
66 
70 
82 
90 
95 

99 

100 
106 


IV 


4.4 Mapping of Explicit to Implicit Representation 
for Quadric Surfaces 

CHAPTER 5: EXPERIMENTAL RESULTS 

5.1 Introduction 

5.2 Median Filtering on Range Images 

5.3 Application of the Recognition Process to the 

Processed Image Data 

5.4 Application of the Rotation Alignment Algorithm 

5.5 Application of Three-Dimensional Discriminant 

Technique 

CHAPTER 6: CONCLUSIONS 

6.1 Overview 

6.2 Advantages of the Recognition Scheme 

6.3 Future Goals and Research Directions 

REFERENCES 

APPENDIX A 

APPENDIX B 

APPENDIX C 

APPENDIX D 

APPENDIX 


107 

112 

112 

113 


156 

162 


166 

171 

171 

174 

175 

176 
180 
187 
190 
193 
198 


v 


LIST OF TABLES 


Table 2-1 Surface classification using the three-dimensional discriminant 

approach 

Table 4-1 Intersection of ellipsoid with planes 

Table 4-2 Intersection of quadric cylinder with planes 

Table 4-3 Intersection of sphere with planes 

Table 4-4 Intersection of quadric cone with planes 

Table 4-5 Intersection of hyperboloid of one sheet with planes 

Table 4-6 Intersection of hyperboloid of two sheets with planes 

Table 4-7 Intersection of elliptic paraboloid with planes 

Table 4-8 Intersection of hyperbolic paraboloid with planes 

Table 4-9 Intersection of the hyperbolic cylinder with planes 

Table 4-10 Intersection of the parabolic cylinder with planes 


Table 4-11 Intersection of the parallelepiped with planes 

Table 4-12 Various curves intercepted by quadric surfaces when 

intersected with the planes z = k and y =k 

Table 4-13 Feature vector (representing the prescence or absence of curves) 

for each of the quadric surfaces 

Table 5-1 Comparison of coefficients evaluated for the original and 

the processed images of a sphere 


18 

55 

62 

64 

73 

82 

88 

95 

96 
100 
103 
106 

108 

109 

159 


vi 


Table 5-2 


Table 5-3 


Table 5-4 


Table 5-5 


Table 5-6 


Table 5-7 


Table 5-8 


Table 5-9 


Table 5-10 


Table 5-11 


Table 5-12 


Table 5-13 


Comparison of coefficients evaluated for the original and 

the processed images of a cylinder 

Curves intercepted by the two planes, z - k and y = k, 
with real raw and processed range data of the sphere 

Curves intercepted by the two planes, z = k and y = k, 
with real raw and processed range data of the cylinder 

Curves intercepted by the two planes, z = k and y = k, 
with real raw and processed range data of the quadric cone 

New coefficients of the raw image data of sphere after 

alignment 

New coefficients of the 3 x 3 median filtered image data of 

sphere after alignment 

New coefficients of the 5 x 5 median filtered image data of 

sphere after alignment 

New coefficients of the raw image data of cylinder 

after alignment 

New coefficients of the 3 x 3 median filtered image 

data of cylinder after alignment 

New coefficients of the 5 x 5 median filtered image data 

of cylinder after alignment 

New coefficients of an unknown simulated data obtained 

after alignment 

New coefficients of an unknown simulated data obtained 
after alignment 


160 


161 


161 


162 


163 


163 


164 


164 


165 


165 


167 


Vll 


Table 5-14 New coefficients of an unknown simulated data obtained 

after alignment 

Table 5-15 Surface characterization using 3-D discriminant approach 
for simulated data 


168 

169 


viii 


LIST OF FIGURES 


Figure 2-1 

Figure 2-2 

Figure 2-3 
Figure 2-4 
Figure 2-5 
Figure 3-1 

Figure 3-2 

Figure 3-3 

Figure 3-4 
Figure 3-5 
Figure 3-6 
Figure 3-7 
Figure 3-8 
Figure 4-1 
Figure 4-2 


Shape of a surface in the vicinity of an elliptic, hyperbolic, 

and parabolic point 

A set of eight view-independent surface types for a visible 

surface 

Two right-handed rectangular coordinate systems 

Relation between the coordinates of P upon translation 

Two rectangular coordinate systems having the same origin . 

Quadric representations of Real ellipsoid. Hyperboloid of 
one sheet. Hyperboloid of two sheets, and Real quadric cone 

Quadric representations of Elliptic paraboloid. Hyperbolic 
paraboloid, Elliptic cylinder, and Parabolic cylinder 

Quadric representations of Hyperbolic cylinder and 

Parallelepiped 

Raw range image of the sphere 

3x3 median filtered image of the raw sphere 

5x5 median filtered image of the raw sphere 

7x7 median filtered image of the raw sphere 

Rotation transformation of the coordinate system 

List of all quadric surfaces 

Ellipsoid detailed view: horizontal intersections 


10 

11 

14 

14 

15 

20 

21 

22 

25 

25 

27 

27 

36 

50 

53 


IX 


Figure 4-3 
Figure 4-4 
Figure 4-5 
Figure 4-6 
Figure 4-7 
Figure 4-8 
Figure 4-9 
Figure 4-10 
Figure 4-11 

Figure 4-12 

Figure 4-13 

Figure 4-14 
Figure 4-15 

Figure 4-16 

Figure 4-17 
Figure 4-18 

Figure 4-19 
Figure 4-20 


Ellipsoid detailed view: vertical intersections 

Cylinder detailed view: horizontal intersections 

Cylinder detailed view: vertical intersections 

Cylinder: lateral view 

Intersection of a plane and a sphere 

Quadric cone detailed view: horizontal intersections 

Quadric cone detailed view: vertical intersections 

Quadric cone: Lateral view 

Hyperboloid of one sheet detailed view: horizontal 

intersections and z = 0 

Hyperboloid of one sheet detailed view: horizontal 

intersections and z = 

Hyperboloid of one sheet detailed view: vertical 

intersections 

Hyperboloid of one sheet: lateral view 

Hyperboloid of two sheets detailed view: horizontal 

intersections 

Hyperboloid of two sheets detailed view: vertical 

intersections 

Hyperboloid of two sheets: lateral view 

Elliptic paraboloid detailed view: horizontal 

intersections 

Elliptic paraboloid detailed view: vertical intersections 
Hyperbolic paraboloid detailed view: horizontal 
intersections 


54 

58 

60 

61 

65 

69 

71 

72 

76 

78 

80 

81 

85 

87 

89 

92 

94 

97 


x 


Figure 4-21 
Figure 4-22 

Figure 4-23 
Figure 4-24 
Figure 4-25 
Figure 5-1 
Figure 5-2 
Figure 5-3 
Figure 5-4 
Figure 5-5 
Figure 5-6a 

Figure 5-6b 

Figure 5-6c 

Figure 5-6d 

Figure 5 -7 a 

Figure 5-7b 

Figure 5-7c 


Hyperbolic paraboloid detailed view: vertical intersections 

Hyperbolic cylinder detailed view: horizontal 

intersections 

Hyperbolic cylinder detailed view: vertical intersections 

Parabolic cylinder detailed view: horizontal intersections 

Parabolic cylinder detailed view: vertical intersections 

Raw range image of the sphere with its background 

Raw range image of the sphere after segmentation 

3x3 median filtered image of the raw sphere 

5x5 median filtered image of the raw sphere 

7x7 median filtered image of the raw sphere 

First derivative with respect to x-axis of the sphere raw 

image 

First derivative with respect to y-axis of the sphere raw 

image 

Second derivative with respect to x-axis of the sphere raw 

image 

Second derivative with respect to y-axis of the sphere raw 

image 

First derivative with respect to x-axis of the sphere filtered 

with a mask size of 3 x 3 

First derivative with respect to y-axis of the sphere filtered 

sphere filtered with a mask size of 3 x 3 

Second derivative with respect to x-axis of the sphere filtered 
with a mask size of 3 x 3 


98 


101 

102 

104 

105 
114 

114 

115 

116 
116 


117 


118 


119 


120 


121 


122 

123 


xi 


124 


Figure 5-7d 

Figure 5-8a 

Figure 5-8b 

Figure 5- 8c 

Figure 5-8d 

Figure 5 -9a 

Figure 5-9b 

Figure 5-9c 

Figure 5-9d 

Figure 5-10 
Figure 5-11 
Figure 5-12 
Figure 5-13 
Figure 5-14 
Figure 5-15 
Figure 5-16 


Second derivative with respect to y-axis of the sphere filtered 

with a mask size of 3 x 3 

First derivative with respect to x-axis of the sphere filtered 

with a mask size of 5 x 5 

First derivative with respect to y-axis of the sphere filtered 

with a mask size of 5 x 5 

Second derivative with respect to x-axis of the sphere filtered 

with a mask size of 5 x 5 

Second derivative with respect to y-axis of the sphere filtered 

with a mask size of 5 x 5 

First derivative with respect to x-axis of the sphere filtered 

with a mask size of 7 x 7 

First derivative with respect to y-axis of the sphere filtered 

with a mask size of 7 x 7 

Second derivative with respect to x-axis of the sphere filtered 

filtered with a mask size of 7 x 7 

Second derivative with respect to y-axis of the sphere filtered 

filtered with a mask size of 7 x 7 

Sign plot for the original range image of the sphere 

Sign plot for the 3 x 3 filtered image of the sphere 

Sign plot for the 5 x 5 filtered image of the sphere 

Sign plot for the 7 x 7 filtered image of the sphere 

Raw range image of the cylinder with its background 

Raw range image of the cylinder after segmentation 

3x3 median filtered image of the raw cylinder 


125 


126 


127 


128 


129 


130 


131 


132 

134 

135 

136 

137 
139 

139 

140 


Xll 


140 


Figure 5-17 5x5 median filtered image of the raw cylinder 

Figure 5- 18a First derivative with respect to x-axis of the cylinder raw 

image 

Figure 5- 18b First derivative with respect to y-axis of the cylinder raw 

image 

Figure 5- 18c Second derivative with respect to x-axis of the cylinder raw 

image 

Figure 5-18d Second derivative with respect to y-axis of the cylinder raw 

image 

Figure 5- 19a First derivative with respect to x-axis of the cylinder filtered 

with a mask size of 3 x 3 

Figure 5- 19b First derivative with respect to y-axis of the cylinder filtered 

with a mask size of 3 x 3 

Figure 5- 19c Second derivative with respect to x-axis of the cylinder filtered 

with a mask size of 3 x 3 

Figure 5-19d Second derivative with respect to y-axis of the cylinder filtered 

with a mask size of 3 x 3 

Figure 5-20a First derivative with respect to x-axis of the cylinder filtered 

with a mask size of 5 x 5 

Figure 5-20b First derivative with respect to y-axis of the cylinder filtered 

with a mask size of 5 x 5 

Figure 5-20 Second derivative with respect to x-axis of the cylinder filtered 

with a mask size of 5 x 5 

Figure 5-20d Second derivative with respect to y-axis of the cylinder filtered 
with a mask size of 5 x 5 


141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 


xm 


153 


Figure 5-21 
Figure 5-22 
Figure 5-23 
Figure 5-24 
Figure 5-25 
Figure 6-1 

Figure 6-2 

Figure 6-3 


Sign plot for the original cylinder 

Sign plot for the 3 x 3 filtered cylinder image 

Sign plot for the 5 x 5 filtered cylinder image 

Best-fit plot for the sphere raw image 

Best-fit plot for the cylinder 

Delta rocket composed of cylindrical and conical 

shapes 

Conical domes and cylinderlike body make up the 

space probe 

Cylinderical space station with a half sphere dome top 


154 

155 

157 

158 


172 


173 

173 


xiv 


CHAPTER ONE 


INTRODUCTION 


1.1 Introduction 

One of the most important tasks in computer vision is that of three-dimensional 
object recognition. Success has been limited to the recognition of symmetric objects. 
Recently, research has concentrated on the recognition of small numbers of asymmetric 
objects as well as objects placed in complex scenes. Unlike the recognition procedure 
developed for intensity-based images, the recent development of active and passive 
sensors extracting quality range information has led to the involvement of explicit 
geometric representations of the objects for the recognition schemes [1, 2]. Location 
and description of three-dimensional objects from natural light images are often 
difficult to determine. However, range images give a more detailed and direct 
geometric description of the shape of the three-dimensional object. A brief introduc- 
tion to range images and the laser range-finder is presented in Section 1.2. In Section 
1.3, a precise global definition of the object recognition problem is discussed. The 
objective of this dissertation and its relevance to the global three-dimensional problem 
is presented in Section 1.4. 

1.2 Range Image and Data Acquisition 

Range images share the same format as intensity images, i.e., both of these 
images are two-dimensional arrays of numbers, the only difference being that the 
numbers in the range images represent the distances between a sensor focal plane to 
points in space. The laser range-finder or tracker [3] is currently the most widely used 
sensor. The laser range-finder makes use of a laser beam which scans the surfaces in 
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the scene of observation from left to right and top to bottom. Thus the distances 
obtained measure both depth and scanning angle. The principle of triangulation is util- 
ized to obtain the three-dimensional coordinate of each pixel. Unless a specific algo- 
rithm demands a special form of the range images, it is usually this depth information 
which is utilized for the recognition process. Active triangulation techniques use an 
extra source of light to project some pattern onto the objects to be measured, thereby 
reducing complexity of the stereo matching problem [4, 5]. Many industrial and navi- 
gational robotic tasks such as target identification and tracking, automated assembly, 
bin picking, mobile robots, etc., will be better accomplished if such explicit depth 
information can be efficiently obtained and accurately interpreted. 

Modeling human vision is a complex process. To date, machine vision systems 
can hardly perform a fraction of the capabilities of the human visual system. An 
efficient mechanism which can acquire relevant information from the three-dimensional 
world and subsequently form models of the real world will, to some extent, bridge the 
gap between machine and human capabilities. 

1.3 Definition of the Object Recognition Problem 

Three-dimensional object recognition is vast problem. In the course of the 
succeeding text, we will give a somewhat precise definition of this problem. 

In the real world, the things human see and feel are primarily solid objects. 
When people view objects for the first time, they attempt to collect information from 
various aspects of the object. This process of collecting and forming information 
about unknown objects is known as model formation [8]. After gaining familiarity 
with many objects, we are able to identify objects from any arbitrary viewpoint 

without further investigation. 

The human vision system has the capability of analyzing and determining 
only the color but also the spatial orientation of objects relative to a fixed coordinate 
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system. Since we are interested in an automatic, computerized process to recognize 
objects, the input data we use must be compatible with available digital computers. 
Hence, two-dimensional matrices of numerical values usually known as digitized sen- 
sor data, constitute the information that is processed to describe or recognize three- 
dimensional objects. The sensor used for this process can be a passive sensor, like a 
camera, or an active sensor, such as a laser range mapper. Summarizing, the three- 
dimensional recognition problem constitutes a detailed completion of model formation 
of the object leading to an in-depth knowledge of its shape and orientation with respect 
to a fixed view of the real world. 

1.4 Objectives and Organization of the Report 

An approach based on two-dimensional analytic geometry to recognize a series of 
three-dimensional objects is presented in this research. Among the various three- 
dimensional objects considered are the hyperboloids of one and two sheets, ellipsoids, 
spheres, circular and elliptical quadric cones, circular and elliptical cylinders, parabolic 
and hyperbolic cylinders, elliptic and hyperbolic paraboloids, and parallelepipeds. 

The difficulties in recognizing three-dimensional objects stems from the complex- 
ity of the scene, the number of objects in the database and the lack of a priori infor- 
mation about the scene. Techniques vary based upon the difficulty of the recognition 
problem. In our case we attempt to recognize segmented objects in range images. 

Location and orientation of three-dimensional objects has always been the most 
complex issue in many computer vision applications. Algorithms for a robust three- 
dimensional recognition system must be view-independent. Herein, we have developed 
a technique to determine the three-dimensional object location and orientation in range 
images. Once the object lies in a desired stable rest position, our proposed recognition 
scheme quickly and accurately classifies it as one of the objects mentioned above. In 
comparison to most of the present day methods utilized for range image object 
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recognition, our proposed approach attacks the problem in a different manner and is 
computationally inexpensive. 

Chapter Two reviews some of the earlier and current work in this area. It 
includes a review of some of the mathematical concepts associated with three- 
dimensional object recognition. A mathematical quadric classification method based 
on a three-dimensional discriminant is discussed while in this chapter. In chapters 
Three and Four we discuss, in detail, our proposed three-dimensional approach. 
Chapter Three addresses the various pre-processings steps involved prior to the appli- 
cation of the recognition algorithm. Median filtering, segmentation, three-dimensional 
coefficient evaluation, and rotation alignment being some of them. The demerits of 
existing schemes in the area of three-dimensional object recognition and the unique- 
ness and improvizations brought about through our recognition procedures are also dis- 
cussed in Chapter Three. In Chapter Four, after a brief discussion of the practical 
merits of using planar intersections, characteristic feature vectors are obtained for each 
of the quadric surfaces under investigation. Results are summarized in Chapter Five. 
A large set of real range images of spheres, cylinders, and cones were utilized to test 
the proposed recognition scheme. Results obtained for simulated data of other quadric 
surfaces, namely, hyperboloids and paraboloids are also tabulated in Chapter Five. 
Chapter Six concludes with a discussion of possible areas for future investigation. 



CHAPTER TWO 


BACKGROUND 


2.1 Introduction 

Past and present research in the field of three-dimensional object recognition is 
reviewed in Section 2.2. Surface curvatures which are widely utilized in this research 
area are briefly reviewed in Section 2.3. Section 2.4 investigates a three-dimensional 
approach to classification and reduction of quadrics as presented by Olmstead [24], 
wherein various invariant features of the quadratic form under translation and rotauon 

are discussed. 

2.2 Literature Review 

Many of the currently available techniques for describing and recognizing three- 
dimensional objects are based on the principle of segmentation. Segmentation is the 
process in which range data is divided into smaller regions (mostly squares) (4). 
These small regions are approximated as planar surfaces or curved surfaces based upon 
the surface mean and Gausssian curvatures. Regions sharing similar curvatures are 
subsequently merged. This process is known as region growing. Other approaches 
[6-10] characterize the surface shapes while dealing with the three-dimensional recog- 
nition problem. Levine et al. [11] briefly review various works in the field of segmen- 
tation, where segmentation has been classified inlo region-based and edge-based 
approaches. Again surface curvatures play an important role for characterization in 

each of these approaches. 

Grimson et al. [12] discuss a scheme utilizing local measurements of three- 
dimensional positions and surface normals to identify and locate objects from a known 
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set. Objects are modeled as polyhedra with a set number of degrees of freedom with 
respect to the sensors. The authors claim a low computational cost for their algorithm. 
Although they have limited the experiments to one model, i.e., data obtained from one 
object, they claim that the algorithm can be used for multiple object models. Also, 
only polyhedral objects with a sufficient number of planar surfaces can be used in their 

scheme. 

Another paper by Faugeras et al. [13] describes surfaces by curves and patches 
which are further represented using linear parameters such as points, lines and planes. 
Their algorithm initially reconstructs objects from range data and consequently utilizes 
certain constraints of rigidity to recognize objects while positioning. They arrive at the 
conclusion that for an object to be recognized, at least a certain area of the object 
should be visible (approx. 50%). They claim their approach could be used for images 
obtained using ultrasound, stereo, and tactile sensors. 

Hu and Stockman [14] have employed structured light as a technique for three- 
dimensional surface recognition. The objects are illuminated using a controlled light 
source of a regular pattern, thereby creating artificial features on the surfaces which are 
consequently extracted. They claim to have solved the problem known as "grid line 
identification." From the general constraints, a set of geometric and topological rules 
are obtained which are effectively utilized in the computation of grid labels which are 
further used for finding three-dimensional surface solutions. Their results infer that 
consistent surface solutions are obtained very fast with good accuracy using a single 
image. 

Recognition of polyhedral objects involves the projection of several invariant 
features of three-dimensional bodies onto two-dimensional planes [15]. Recently, 
recognition of three-dimensional objects based upon their representation as a linear 
combination of two-dimensional images has been investigated [16]. Transformations 
such as rotation and translation have been considered for three-dimensional objects in 
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terms of the linear combination of a series of two-dimensional views of the objects. 
Instead of using transformations in three-dimensions, it has been shown that the pro- 
cess is the equivalent of obtaining two-dimensional transformations of several two- 
dimensional images of the objects and combining them together to obtain the three- 
dimensional transformation. This procedure appears computationally intensive. 

Most of the techniques and algorithms mentioned above have a common criterion 
for classifying the three-dimensional objects in the final phase. They have a database 
of all the objects they are trying to recognize and hence try to match features from the 
test samples to the features of the objects in the database. 

Fan et al. [17] use graph theory for decomposing segmentations into subgroups 
corresponding to different objects. Matching of the test objects with the objects in the 
database is performed in three steps: the screener, which makes an initial guess for 
each object; the graph matcher, which conducts an exhaustive comparison between 
potential matching graphs and computes three-dimensional transformation between 
them; and finally, the analyzer, which based upon the results from the earlier two 
modules conducts a split and merge of the object graphs. The distinguishing aspect of 
this scheme is that the authors used occluded objects for describing their proposed 
method. 

As has been mentioned, most of the present research on three-dimensional objects 
utilize range imagery rather than stereo images. But at the same time, it should be 
noted that it was stereo imagery which, to a large extent, was initially used to investi- 
gate the problem of three-dimensional object recognition. 

Forsyth et al. [18] use stereo images to obtain a range of invariant descriptors in 
three-dimensional model-based vision. Initially, they demonstrate a model-based 
vision system that recognizes curved plane objects irrespective of the pose. Based 
upon image data, models are constructed for each object and the pose is computed. 
However, they mainly describe three-dimensional objects with planar faces. 
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Lee and Hahn [19] have actually dealt with an optimal sensing strategy. Their 
main objective is to obtain valuable and effective data or information from three- 
dimensional objects, which subsequently could be used to describe and recognize 
natural quadric surfaces. Other works on stereo vision can be found in references 20, 

21, 22 and 23. 

The visible-invariant surface characteristics mentioned before are the Gaussian 
curvature (K) and the mean curvature (H), which are referred to collectively as surface 
curvatures. Mean curvature is an extrinsic surface property, whereas Gaussian curva- 
ture is intrinsic. In the following section we briefly describe these two widely used 
invariant surface characteristics for three-dimensional objects. 

2.3 Differential Geometry of Surfaces: Mean and Gaussian Curvatures 

Mean and Gaussian curvatures [8] are identified as the local second-order surface 
characteristics that possess several desirable invariance properties and represent extrin- 
sic and intrinsic surface geometry, respectively. The explicit parametric form of a gen- 
eral surface S in E 3 (three-dimensional Euclidean space) with respect to a known 

coordinate system is given as 

S = j(x(u,v), y(u,v), z(u,v)): (u,v) e d|, ( 21 ) 

where D is any surface patch and is a subset of E . 

However if the depth maps are assumed to be in the form of a graph surface 
(Monge patch surface) [8], then S can be written as 

S =|(x,y,z(x,y)), (x,y) £ d|, 

where z(x,y) is the depth at a point (x,y) in a given range image. 
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The representations for the Gaussian and the mean curvatures are as follows: 


Gaussian curvature, K, is defined by, 

9 1 2 3 4 z 3 2 z 

3x 2 dy 2 


d 2 z 


dxdy 


l + (^) 2 + (^) 2 


dx 


dy 


( 2 . 2 ) 


Mean curvature, H, is defined by, 


(2.3) 


d 2 z ! 8 2 z ( 8 2 z 

3zl 

d 2 z 

r 8y 2 

r 8zl 

dx 

i J 

3z dz d 2 z 

'dx'dy dxdy 

8x 2 dy 2 dx 2 

[8y^ 

H 

F 

dz 

dx 

L J 

+ 

«■ 

dz 

dy 

^3/2 


Both of these curvatures are invariant to translation and rotation of the object as long 
as the object surface is visible. 

Based upon the sign of the Gaussian curvature, individual points in the surface 
are locally classified into three surface types as shown in Figure 2-1: 

(a) K > 0 implies an elliptic surface point, 

(b) K < 0 implies a hyperbolic surface point, and 

(c) K = 0 implies a parabolic surface point. 

Besl and Jain in their paper [8] have shown that the Gausssian and mean curva- 
tures together can be utilized to give a set of eight different surfaces as shown in Fig- 
ure 2-2: 


1) H < 0 and K > 0 

2) H > 0 and K > 0 

3) H < 0 and K = 0 

4) H > 0 and K = 0 


implies a peak surface, 
implies a pit surface, 
implies a ridge surface, 
implies a valley surface 
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(b) Hyperbolic point (K < 0) 



(a) Elliptic point (K > 0) 


(c) Parabolic point (k = 0) 


Figure 2-1. Shape of a surface in the vicinity of an elliptic, hyperbolic, and parabolic 
point. 



Ridge Surface H < 0, K = 0 


Valley Surface H > 0, K = 0 


Figure 2-2. A set of eight view-independent surface types for a visible surface 
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5) H = 0 and K = 0 

6) H = 0 and K < 0 

7) H < 0 and K < 0 

8) H > 0 and K < 0 


implies a flat surface, 
implies a minimal surface, 
implies a saddle ridge surface, 
implies a saddle valley surface. 


2.4 Three-Dimensional Discriminant 

In this section we investigate a three-dimensional approach to classification and 
reduction of quadrics as presented by Olmstead [24], which looks into the invariants of 
the quadratic form under translation and rotation of three-dimensional objects. 

The general quadric surface of second degree in the three variables x, y, and z 
can be written in the form 


F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d - 0 
Associated with F(x,y,z) are two matrices: e and E, where 


and 


e = 


a h g 
h b f 
£ f C 


Ip q r d] 

Let the determinant of E be denoted by A, and the determinant of e be denoted 
by D. Also let the cofactors of each element of A be denoted by the corresponding 
capital letters. Three-dimensional surfaces are classified as singular or non-singular, 
based upon E being singular or non-singular. Examples of non-singular surfaces are 
ellipsoids, hyperboloids, and paraboloids. The other quadrics are singular. 
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Let us now consider the two basic transformations, namely translation and rota- 
tion, and try to arrive at some invariant features. Consider the two rectangular right- 
handed coordinate systems as shown in Figure 2-3. Any point in space has two sets of 
coordinates, one for each set of axes. The problem is to find a relationship between 
these two sets of coordinates so that one can convert from one coordinate system to 

the other. 


2.4.1 Translation 

Inspecting Figure 2-4, we see that the coordinates of O' and P in the xyz sys- 
tem are (x 0 ,y o ,z 0 ) and (x,y,z), respectively, and the coordinates of P in the x'y'z' sys- 
tem are (x',y',z'). The two sets of coordinates of P are related by the following trans- 
lation equations: 


x = x' + x Q . 

y = y' + y 0 - 

z = z' + z 0 . 


(2.4) 

(2.5) 

(2.6) 


The set of equations, (2.4), (2.5), and (2.6) relate the coordinates of a point in the 
x'y'z' system to its coordinates in the xyz system. Direct substitution of equations 
(2.4) - (2.6) into F(x,y,z) results in the following theorem: 


Theorem 2.1. For any quadric surface, the coefficients of the second degree terms 
and therefore the matrix e, are invariant under translation. 


2.4.2 Rotation 

Consider the two rectangular coordinate systems as shown in Figure 2-5. With 
respect to the x'y'z' system, let the direction cosines of the x, y, and z axes be 
(Ml.Vj), (X 2 .D 2 .V 2 ). and (X 3 .D 3 .V 3 ), respectively. Then with respect to the xyz sys- 
tem, the direction cosines of the x', y', and z' axes are (Xj, X 2 , X 3 ), (t>i, *> 2 , t> 3 ), and 
(v lt v 2 , v 3 ), respectively. 
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For any point, P, whose coordinates in the two systems are (x,y,z) and (x ,y ,z ), 
the following two sets of rotation equations are obtained: 

x = Xjx' + tijy' + Vjz', 
y = X^x + t) 2 y' + V 2 Z > 
z = A, 3 x' + t>3y' + V3Z', 

and 


x = A,]X + X 2 y + X^z, 
y' = t)jX + o 2 y + O 3 Z, 
z = VjX + v 2 y + V 3 Z, 

which gives rise to the rotation matrix 


A = 


X* v i 
X 2 v 2 v 2 , 
h u 3 v 3 


(2.7) 


where the elements of the rows (or columns) are direction cosines of perpendicular 
directions. Direct calculation results in the following theorem: 


Theorem 2.2. The determinant D of the rotation matrix A is equal to 1. 


Before arriving at a particular set of invariant features of a quadric, we first 
describe a plane of symmetry of a certain type, called a principal plane. 

Definition 2.1 A principal plane is a diametrical plane that is perpendicular to the 

chord it bisects [24]. 

Consider the matrix e again: 


e = 


a h g 


h b f . 

Ifi f cj 


The eigen-values of the matrix e can be calculated from 
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a - k h 

h b - k 

g f 


g 

f 

c - k 


= 0 . 


This cubic equation in k is called the characteristic equation of the matrix e. Its 
roots are called the characteristic roots of e. The quantities given below are found to 
be invariant as a consequence of the following theorem [25]. 


Theorem 2.3 If the second degree equation F(x,y,z)=0 is transformed by means of a 
translation or a rotation with fixed origin, the following quantities are invariant: 

D, A, p 3 , p 4 , I, J, kj, k 2 , and k 3 , where D, A are the determinants of the matrices e 

and E, respectively; and p 3 and p 4 are the ranks of the matrices e and E, respec- 
tively. Also 


I = a + b + c. 


J = ab + ac + be - f 2 - g 2 - h 2 , 

and finally kj, k 2 , and k 3 are the characteristic roots of e. 

Based upon the above set of invariants, surface classifications are listed in Table 

2 - 1 . 

In Chapters Three and Four, we discuss our proposed recognition scheme in 


detail. 
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Number 

Surface 

P3 

P4 

Sign of A 

k’s same sign 

1 

Real ellipsoid 

3 

4 

- 

yes 

2 

Hyperboloid of one sheet 

3 

4 

+ 

no 

3 

Hyperboloid of two sheets 

3 

4 

* 

no 

4 

Real quadric cone 

3 

3 


no 

5 

Elliptic paraboloid 

2 

4 

~ 

yes 

6 

Hyperbolic paraboloid 

2 

4 

+ 

no 

7 

Real elliptic cylinder 

2 

3 


yes 

8 

Hyperbolic cylinder 

2 

3 


no 

9 

Parabolic cylinder 

1 

3 




Table 2-1. Surface classification using the three-dimensional discriminant approach. 
p 3 is the rank of matrix e and p 4 is the rank of matrix E. The characteristic roots of 
the matrix e are referred by k’s. 




CHAPTER THREE 


QUADRIC SURFACE REPRESENTATION 


3.1 Introduction 

Section 3.2 considers the various three-dimensional quadric surfaces used in the 
recognition process. While describing each of these objects, we will be considenng 
the surfaces with their centers aligned to the origin of our coordinate system. Section 
3.3 explains our quadric recognition algorithm in detail. This section also addresses 
the acquisition of range data and the necessary pre-processing steps, the representation 
of quadric surfaces by a second degree polynomial, and the rotation alignment algo 
rithm whereby each of the quadric surfaces are placed in a coordinate system of our 
choice. The merits of the proposed technique are addressed while considenng the 
improvisations brought about in the recognition of three-dimensional objects (espe- 
cially quadrics) in Section 3.4. 


3.2 Quadric Surface Description 

In this section by means of Figures 3-1, 3-2, and 3-3, we illustrate and represent 
the following three-dimensional quadric surfaces which are considered for the recogni- 
tion process: ellipsoids, the hyperboloids of one and two sheets, quadric cones, elliptic 
paraboloids, hyperbolic paraboloids, elliptic cylinders, hyperbolic cylinders, parabolic 

cylinders, and parallelepipeds. 


Most three-dimensional objects of practical use consist of at least one of the sur- 
faces described above. All the representations of surfaces which were described above 
hold true under ideal conditions, he., when the source daia is perfect, exact pose and 


orientation of the objects are known, the system is noiseless, etc 


However in the real 
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X 2 V 2 z 

Real Ellipsoid: — + — + — 

3. D C 


Hyperboloid of one sheet. — 2 + c2 



Figure 3-1. Quadric representations of Real ellipsoid. Hyperboloid of one sheet 
Hyperboloid of two sheets, and real quadric cone. 




Elliptic cylinder: JL + -^ = 1 Parabolic cylinder, x 2 + 2rz - 0 

Figure 3-2. Quadric representations of Elliptic paraboloid, Hyperbolic paraboloid 
Elliptic cylinder, and Parabolic cylinder. 
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world, practically none of these conditions hold true. Any set of data, whether it is 
derived or generated from a passive (camera) or an active sensor (laser range mapper), 
can at best be approximated to a second-degree polynomial. Whether this polynomial 
accurately represents a surface or not, and if so, how these coefficients (representa- 
tion) can be chosen to come close to recognizing a three-dimensional object, is the 
whole issue of the recognition problem. 

In the next few sections, while formulating our recognition scheme, we describe 
one such technique which generates ten coefficients (which are sufficient under ideal 
conditions) to describe all the objects of interest [26]. 

However, before elaborating on the recognition scheme, an overview of the tech- 
nique is presented. The recognition scheme utilizes a two-dimensional discriminant 
(which is a measure for distinguishing two-dimensional curves) to recognize three- 
dimensional surfaces. Instead of utilizing the ten generated coefficients and attempting 
to recognize the surface from its quadric representation, the quadrics are identified 
using the information resulting from the intersection of the surface with different 
planes. If the surface is one of those listed above, there are five possible two- 
dimensional curves that may result from such intersections, (i) a circle, (ii) an ellipse, 
(iii) a parabola, (iv) a hyperbola, and (v) a line. Thus, a feature or pattern vector 
with five independent components can be formed for characterizing each of the sur- 
faces. 

3.3 Recognition Scheme 

Our recognition scheme consists of the following steps: 

(1) acquisition of the range data and conducting the pre-processing steps, 

(2) description and representation of objects as general second degree surfaces, 

(3) determination of the location and orientation of the objects with respect to a 
desired coordinate system, 
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(4) performance of the rotation and translation transformations of the object so as to 

place it in a stable desired coordinate system, 

(5) use of the principle of two-dimensional discriminants to classify the various curves 

obtained by intersecting the surfaces with planes, and 

(6) acquisition of an optimal set of planes sufficient enough to distinguish and 
recognize each of the quadric surfaces. Angular bounds within which every 
surface yields a distinct set of curves are determined in step 6. 

The range data, as mentioned in Chapter One, is a pixel-by-pixel depth 
from the point of origin of the laser to the point where the beam impinges on a sur- 
face. The objects are scanned from left-to-right and top-to-bottom. A grid frame may 
consist of 256 x 256 pixels. Before this range data is applied to the object classifier, it 

has to undergo the following pre-processing steps. 

a) median filtering, and 

b) segmentation. 

3.3.1 Median Filtering 

Conventionally, a rectangular window of size M x N is used in two dimensional 
median filtering. As in our case [27], experiments were performed with square win- 
dows of mask sizes 3 x 3 and 5 x 5. Salt and pepper noise in the range images used 
in this research was uniformly distributed throughout. Irrespective of the mask size, 
the range information at every pixel in the image is replaced by the median of the pix- 
els contained in the M x M window centered at that point. Referring to Figure 3-4 
and keeping in mind that the black pixels correspond to the background and the white 
pixels to the object, black pixels inside the object are referred to as pepper noise and 
white pixels in the black background are referred to as salt noise. Figure 3-5 is 
obtained as a result of a 3 x 3 mask being moved over the entire image. The image 
looks as sharp as the original image though some of the noise still exists. A 5 x 5 and 



Figure 3-4. Raw range image of the sphere. 



Figure 3-5. 3x3 median filtered image of the raw sphere. 
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a 7 x 7 mask removes all of the visual salt and pepper noise, but the images as seen in 
Figures 3-6 and 3-7 respectively, to some extent, have lower contrast than the original 

image. 

Once a range image is filtered using a median filter of different masks, the next 
concern is to study the changes to the original data which have been brought about by 
filtering. Evaluating curvatures is one good way of distinguishing similarities and dis- 
similarities among the filtered images and the original range data. 

First and second order derivatives are evaluated along the x and y axes to check 
the uniformity of the original and the filtered images. Approximating, the first-order 
derivative for a pixel (Ay) centered at i j is given as: 

= -~[(A} + l,j + l - Ay^i) + (A i+ ij - Ay)] 

and 

■3^- = ^-[(A i+ ij+i - A i+1 j) + (Ay+i “ Ay)], 

dy 2t 

Similarly approximating, the second order derivatives for a pixel centered at Ay is 
given as: 

= -V [A Mj - 2A io + A w.J ] 

9x z £ 

and 

~ ~[Ay-i - 2Ay + Ay+i], 

dy £ 

where £ represents the spacing between picture cell centers. 

A sign map, which shows the relationships among two neighboring pixels with 
respect to the depth value, was also generated to check the effect of median filtering 
on the original data. Sign maps of some of the experimented quadric surfaces are 


illustrated in Chapter Five. 


Figure 3-6. 5x5 median filtered image of the raw sphere. 



Figure 3-7. 7x7 median filtered image of the raw sphere. 
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3.3.2 Segmentation 

Since isolated objects instead of complex scenes are considered, a simple thres- 
holding whereby the object is separated from the background is utilized for the seg- 
mentation process. In the case where objects are irregular or a scene consists of a 
cluster of objects, Gaussian and mean curvatures have to be utilized to sub-divide the 
scene into planar or curved surfaces. Each surface is then recognized separately. 
Range image segmentation has been extensively studied by Levine et al. [9]. 

Now that the available range data has been processed to eliminate salt and pepper 
noise, we can now utilize the image data to obtain the quadric surface which best fits 
the data. To achieve this goal, we need to determine the coefficients of a second 
degree polynomial representation for the three-dimensional surface. 

3.3.3 Three-Dimensional Coefficients Evaluation 

Our objective is to obtain a surface described by Equation (3.1) from a given set 
of data (range) points. We assume that the data is a set of range-image samples 
obtained from a single surface which can be described by a quadric equation. 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0. (3.1) 

We shall therefore define the best description to be the one which minimizes the 
mean-squared error (MSE) between the range data and the quadric [26]. 

Equation (3.1) in vector notation becomes 

F(x,y,z) = a T p = 0, 2 ) 

where a T = [ a b c 2f 2g 2h 2p 2q 2r d ] and p T = [ x 2 y 2 z 2 yz zx xy x y z 1 ]. 

The error measure for any data point (x,y,z) can be measured by evaluating 
F(x,y,z). If this point lies exactly on the surface then, F(x,y,z) = 0, meaning that the 


error is zero. 
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The mean-squared error, E, is defined as 

E = min £llFll 2 . (3.3) 

• s 

In vector notation, Equation (3.3) becomes 

E = min £a T pp T a = min a T Ra, (3 4) 

a s a 

where R is the scatter matrix for the data set equal to 

R = £PP T - (3-5) 

s 

Minimizing E leads to a trivial solution of a = 0, implying all the coefficients are zero. 
We therefore attempt to find the minimum of a^Ra with respect to a, subject to some 
constraint K(a) = k. 

Let 


G(a) = a T Ra (3-6) 

and 

K(a) = a T Ka, (3-7) 

where K is another undetermined constant matrix. Using Lagrange’s method [28], we 
write the function 


U = G(a) - XK(a), (3-8) 

where X again is an undetermined constant. To find a minimum solution for U, we 
form 


au 

da 


= 2(R - XK)a = 0. 


(3.9) 
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Solving = 0 and K(a) = k simultaneously, we find a and X to give a minimum 
9a 

solution. We wish to evaluate the constraint K(a) such that it gives a non-zero solu- 
tion for a for all the quadric surfaces of interest. 

In order to determine the function of the coefficient vector a which is invariant to 
translation and rotation, we write the quadric equation as 

F(x,y,z) = F(v) = v‘Dv + 2v l q + d = 0, (3.10) 

where 


v = 



(3.11) 


D = 


a h g 
h b f 
.g f c 


( 3 . 12 ) 


and 


q = 


p 

q • 

.r 


(3.13) 


After carrying out the transformations, translation and rotation, it is observed that the 
second-order terms and the eigen-values are the only invariants of D under translation 


and rotation, respectively. 


We now derive a function of the eigen-values of D, i.e., f(X), which will allow us 
to obtain all of the quadrics of interest. The constraint should be in a quadratic form, 
such that when we substitute it in 


— = 2(R - MC)a = 0, (3-14) 

9a 


we get a linear equation from which we can solve for a. 



31 


From reference 29, a good choice for the constraint f(X.) is 


fto = Etf = 1. 


(3.15) 


i.e., 


= tr(D 2) = a 2 + b 2 + c 2 + 2f 2 + 2g 2 + 2h 2 

T 

Writing it in the form of equation K(a) = a Ka: 


tr(D 2 ) = a T 


K 2 0 
0 0 


a. 


where the constraint matrix K 2 is 


K, = 


1 0 0 0 0 0 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1/2 0 0 

0 0 0 0 1/2 0 

0 0 0 0 0 1/2J 


Equation Ra = XKa, can now be written as 


(3.16) 


(3.17) 


(3.18) 


C B 
B t A 


= X 


k 2 0 

0 0 


(3.19) 


where C is the 6 x 6 scatter matrix for the quadratic terms a, b, and c; B is the 6 x 4 
scatter matrix for the mixed terms 2f, 2g, and 2h and A is the 4x4 scatter matrix for 
the linear and constant term, i.e., 2p, 2q, 2r, and d. P is the 6x1 vector of the qua- 
dratic coefficients and a is the 4 x 1 vector of the linear and the constant coefficients. 

Solving Equation (3.19), we get 


and 


Cp + Ba = XK 2 p 


( 3 . 20 ) 


32 


B T p + Aa = 0. 


(3.21) 


From Equation (3.21) we get 


a = - A" 1 B T p. 


( 3 . 22 ) 


Substituting a in Equation (3.20), we have 


( C - BA _1 B t )p = A.K 2 p. 


(3.23) 


Labeling ( C - BA 1 B T ) as M, we have 


Mp = X.K 2 P, 


(3.24) 


which appears similar to an eigen-value problem. Writing 


K 2 as H 2 , where, 


H = 


1 0 0 
0 1 0 
0 0 1 
OP 0 0 
0 0 0 
0 0 0 


0 0 

0 0 

0 0 

1/a/2 0 

0 M<2 

0 0 


0 

0 

0 

0 ’ 
0 

1A/2 


(3.25) 


1 0 0 0 0 0 

0 1 0 0 0 0 

„-i _ 0 0 1 0 0 0 

~ 0 0 0 a/2 0 O' 

0 0 0 o a/2 0 

0 0 0 0 0 a/2 


(3.26) 


We can write MP = A.K 2 P as MP = XHHp, or H 'MH ’HP - XHp. 

Let p' = HP and M' = where M' is a real symmetric matrix, then 


M'p' = XB'. 


(3.27) 


M' has six ^’s and six corresponding Bj’s. 
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For the minimum error solution, we choose the eigen-vector corresponding to the 
smallest eigen-value, i.e., 


Pi = H -1 B'j. (3-28) 

Solving for ctj = -A _1 B T Bi, we have our solution. 

Once the procedure described in Section 3.3.3 has been performed, the median 
filtered range data can be described as 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0, (3.29) 

where the values of the coefficients a, b, c, f, g, h, p, q, r, and d are known. Generally 
speaking, all of the objects in the experiments generate all ten coefficients as is shown 
in Chapter Five. The question now is: How can we distinguish one object from the 
another and how accurately can we describe the recognized object? In the following 
sections of this chapter and Chapter Four, we describe the necessary scheme to solve 
the recognition problem of quadric surfaces. 

3.3.4 Evaluation of the Rotation Matrix 

The determination of the location and orientation of a three-dimensional object is 
one of the central problems in computer vision applications. It is observed that most 
of the methods and techniques which try to solve this problem require considerable 
pre-processing such as detecting edges or junctions, fitting curves or surfaces to seg- 
mented images and computing high order features from the input images. Since 
three-dimensional object recognition depends not only on the shape of the object but 
also the pose and orientation of the object as weU, any definite information about the 
object’s orientation will aid in selecting the right features for the recognition process. 

In this research we suggest a method based on analytic geometry, whereby all the 
rotation parameters of any object placed in any orientation in space are determined and 
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eliminated systematically. With this approach we are in a position to place the three- 
dimensional object in a desired stable position, thereby eliminating the orientation 
problem. We can then utilize the shape information to explicitly represent the three- 
dimensional surface. 

Any quadric surface can be represented by Equation (3.29) in terms of a second 
degree polynomial of variables x, y, and z. 

Let (x,y,z) describe the coordinates of any point in our coordinate system. As 
shown in Figure 3-8(b), consider a rotation of angle a about the z axis, i.e. in the 
xy-plane. Then the new coordinates in terms of the old are represented as 


x = x'cosa + y'sina 


and 


i.e., the rotation matrix is 


y = -x'sina + y'cosa; 


Ra = 


cosa sina 0 
-sina cosa 0 
0 0 1 


Next, as shown in Figure 3-8(c), consider a rotation about the x axis by an angle p, 
i.e., in the y'z plane, of the same point. The resultant coordinates and the old coordi- 
nates are now related by the following equations: 


y' = y"cos(i + z'sinp 
and 


z = -y"sinP + z'cosP, 


where the rotation matrix is 
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1 0 0 
0 cosp sinP . 
0 -sinP cosP 


Finally as shown in Figure 3-8(d), consider a rotation about the y" axis by an angle 
y, i.e., in the xY plane, then 


and 


z = z"cosy + x"siny 


x' = -z"siny + x"cosy. 


The rotation matrix for the above transformation is 


Observing that 


we obtain the following: 


cosy 0 -siny 


r t = 


0 1 
siny 0 


0 

cosy 


X 


X 

y 

— R a RpRy 

// 

y 

[z J 


x = x"(cosacosy + sinasinpsiny) + y"sinacosp + z"(-sinycosa + cosysinasinp), 

y = x"(— cosysina + sinysinpcosa) + y"cosPcosa + z"(sinysina + cosysinpcosa) 
and 

z = x"sinycosP - y"sinP + z"cosycosp. 

After substituting the new x, y, and z coordinates into Equation (3.29), we get an 
entire set of new coefficients for x" 2 , y" 2 , z" 2 , y”z", x"z", x"y", x", y", and z . 
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These new coefficients are listed below. 


a" = cos 2 "/ [a-cos 2 a + b-sin 2 a] + sin 2 Psin 2 y [a-sin 2 a + b-cos 2 a] 

+ 2sinasinpsinycosacosy (a - b) + c-sin^os 2 ^ 

- sin2a [h-sin^sin 2 ^ - sin2y [f-sinacosp + g-cosacosP 
+ h-cos 2 asinP - h-sinpsin 2 a ]+ sin2Psin 2 Y (fcosa - gsina) + h-sin2acos z y. 


b" = (a-sin 2 a + bcos 2 a)cos 2 p + c-sin 2 P + sin2P [-f-cosa g sina ] 
+ h-sin2acos 2 p. 


o" — sin^y (a-cos 2 a + b-sin 2 a) + (a-sin 2 a + b-cos a)cos 2 ysin P 

+ 2sinasinPsin'>cosacosY (a - b) + c-cos 2 ycos 2 P + sin2a [h-cos 2 7 sin 2 p - h-sin^l (3.32) 
+ cos 2 Ysin2p [f-cosa + g sina] + sin2y [-f-sinacosP + g-cosacosp+ h-cos2asinp]. 


2f" = [(b-cos 2 a + a-sin 2 a + h-sin2a - c)sin2P + (2g-sina + 2f-cosa)cos2P 
+ j^((a — b)sin2a + 2h-cos2a)cosP — (2g-cosa — 2f-sina)sinPj siny. 


cosy 


(3.33) 


2g" = sin2y 


-cos 2 a(a + b sin 2 P) - sin 2 a(a-sin 2 P + b) - c-cos 2 P 


- sinPcosp(f-cosa + g-sina) + h-sinacosacosP 
+ 2cospcos 2 y(f-sina — g-cosa) + 2h-sinP(sin 2 asin2y — cos^acos^y). 


(3.34) 


cosacosp(b - a) - h-sinpsinycospj + sin2psiny(a-sin 2 a - b-cos 2 a + c) 


2h" = sin2a 

+ cosysinP(2g-cosa - 2f-sina) + sin 2 Psiny(2g-sina + 2f-cosa) 
- 2h-cos 2 acosy:osp. 


(3.35) 


2p" = 2cosy [-p-cosa + q sina] - 2sinPsiny [p-sina + q-cosa] - 2r sinycosP. (3.36) 
2q" = 2cosp [p-sina + q-cosa] - 2r-sinp. ^ 3 - 37 ^ 

2r" = 2cosysinp [p-sina + q-cosa] + 2siny [p-cosa - q sina] + 2r-cosycosp. (3.38) 


d" = d. 


(3.39) 
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As seen from the above expressions, all of the coefficients are affected by the rotations 
a, p, and y except for the constant d". 

In order to eliminate the product terms 2f", 2g", and 2h , expressions (3-33) 
(3.35) must be set equal to zero and solved simultaneously. As seen from these three 
expressions, each of them is a function of the rotation angles ft, p, and y. It is not 
possible to analytically find the rotation angles which eliminate the product terms. 
Instead, in the next section we present an iterative technique which performs the elimi- 
nation of the product terms. 

3.3.5 Product Terms Elimination Method 

The product terms yz, xz, and xy in F(x,y,z), denote the rotation terms which are 
to be eliminated. Elimination of all these rotation terms will place the three- 
dimensional surface on a coordinate system plane parallel to our coordinate system. 

Observe that in the presence of a single rotation term, say yz, Equation (3.29) 
takes the form 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2px + 2qy + 2rz + d = 0. 

The equation of the trace of the surface in the yz plane is obtained by setting x = 0. 
An appropriate rotation about the origin in the yz plane by an angle p will eliminate 
the yz term. 

However, in the presence of two or more rotation terms, trying to eliminate a 
second rotation term will force the previously eliminated rotation term to reappear. 
Therefore, there will be at least two rotation terms present. The approach we propose 
is an iterative process, whereby at each stage the object is rotated in each of the coor- 
dinate planes, sequentially. The procedure is repeated until all the product terms are 
eliminated, i.e., the coefficients f, g, and h converge to zero in the limit. 
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Since our aim is to eliminate the rotation terms xy, yz, and xz, let’s exclusively 
consider the coefficients of these rotation terms, namely f, g, and h evaluated in Sec- 
tion 3.3.4. In our iterative procedure we are able to eliminate all of the product terms. 
For example, suppose we wish to eliminate the term xy. By a specific rotation of a 
about the z axis, we will be able to accomplish our goal. However, while executing 
this process, the orientation of the object about the two planes yz and zx, i.e., the 
angles the object make with these two planes have been changed. If we wish to elim- 
inate the yz term, the object has to be rotated about the x axis by an angle (3. How- 
ever, in this instance, while performing the process, the previously eliminated xy term 
reappears though the magnitude of its present orientation has been reduced. Hence by 
iterating the above process, an instance occurs when all the coefficients of the product 
terms converge to zero in the limit. 

Consider the Equations (3.33), (3.34), and (3.35). First eliminate the coefficient h, 
i.e, the xy term. This can be accomplished by rotating the object about the z axis by 
an angle a, whereas [3=7=0. Under these circumstances the new coefficients are as 
shown below. 


2f n = 2g-sina 1 + 2fcosa!, 


2g n = 2g-cosa! - 2f sina 1 , 


and 


2h n = (a - b)sin2a, + 2h-cos2aj = 0, 


where cot2a! = 


b - a 
2h 


As seen above, the coefficient h has been forced to 0. The first digit of the subscript 
refers to the iteration number, whereas the second digit of the subscript denotes the 
number of times the object has been rotated by a specific angle. The remaining 
coefficients a, b, c, p, q, and r also reflect changes brought about by the above rotation. 
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The new coefficients are 

a n = a cos 2 a! + bsin 2 ^ - 2hsina 1 cosa 1 , 
b n = b-cos 2 ^ + a-sin 2 ^ + 2hsina 1 cosa 1 , 

c n = c > 

2p n = 2p-cosa! - 2q-sina 1 , 

2cjj| — 2p-sinai 4" 2(j'Cos(X| s 
and 


2f] i - 2r. 

The new quadric equation is 

F(x,y,z) = a n x 2 + b u y 2 + c„z 2 + 2f„yz + 2g u xz + 2 Pll x + 2q„y + 2r„z + d = 0. 

Consider the second step wherein the coefficient corresponding to the yz term is 
forced to zero. In this particular case, the object has to be rotated by an angle p 
about the x axis, where a=Y=0. Under these circumstances, the new rotation 
coefficients (signifying the product terms) become 


2f i2 = (b 12 - c 12 )sin2p! + 2f 11 cos2[} 1 = 0, 


where cot2(5j = 


c H - b 


li 


2f 


ll 


2 Ei2 = 2gii-c°sp 1 , 


and 


2h 12 = ^gn-sinp,. 
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At the same time the other coefficients become 


a 12 _ a ll> 

= Cjj'sin^Pj + b 11 cos 2 p 1 — 2fjj'sinPjCospi, 


c i2 = bu’si^Pj + c 11 -cos 2 Pi + 2f 11 -sinP 1 cospi, 


2pi2 - 2pn, 


2q 12 = 2q 11 -cosp 1 - 2r 11 -sinp 1 , 


and 


2r 12 = 2q 11 sinP 1 + 2r 11 cosp 1 . 

The new quadric equation is: 

F(x,y,z) = a 12 x 2 + b 12 y 2 + c 12 z 2 + 2g 12 xz + 2h 12 xy + 2p 12 x + 2q 12 y + 2r 12 z + d = 0. 

In the final step of the initial iteration, the coefficient corresponding to the xz term is 
forced to zero. In this case, the object is to be rotated by an angle y about the y 
axis, whereas a=p=0. Under these circumstances, the new rotation coefficients beco- 
men 


2fi3 = 2h 12 siny 1 = -^gn-sinpjsiny!. 


2go =(a 13 - c 13 )sin2y 1 + (2g 11 -cosa 1 - fn-sincqjcosp^osly! = 0, 


, . c 12 “ a 12 

where cour{ x = — 

2 El2 


and 


2h 13 = Ih^cosY! = ^gn-sinp^osy!. 
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Let’s now carefully analyze the coefficients of xy, yz, and zx obtained in the final step 
of the first iteration. Consider, for instance, the coefficient corresponding to the yz 
term. It is observed that while proceeding from one step to the other, the new 
coefficients are getting multiplied by the sine or cosine of the concerned angle. This 
implies that in every succeeding step these coefficients are decreasing in their magni- 
tude. To justify the above statement, let us now consider all the coefficients obtained 
in the second iteration. 

At the end of stage 1 of the second iteration, the rotation coefficients become 

2f 21 = 2f 13 -cos02 = -2g 11 -sinP 1 siny ] cos(X2, 

2g 21 = -2f 13 -sina2 = 2g 11 sinP 1 siny 1 sina 2 , 
and 

b 13 - a 13 

2h 2 i = 0, where cot2a 2 = — — • 

zn 13 

At the end of the second stage of the second iteration, the rotation coefficients 
become 


C 2 j — ^21 

2f 22 = 0 where cot2P 2 = — 

2t 2 i 

2g22 = 2g 11 sinp 1 siny 1 sina 2 cosp 2 , 


and 


2h 22 = ^gn-sinPisiny^ino^sinP^ 

Similarly at the end of the final stage of the second iteration, the rotation coefficients 
reduce to 


2f 23 = -2g 11 sinP 1 siny 1 sina 2 sinP 2 siny 2 , 
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2g 23 = 0 


where cot 2 a 2 = 


bp ~ ^13 

2h n ’ 


and 


2h 23 = -2g t 1 sinP 1 siny 1 sina 2 sinP 2 cosY 2 

The terms a 2 , P 2 , and y 2 are the respective rotation angles along the z, x, and y axes 
in the second iteration. Hence it is observed with each iteration that the rotation 
coefficients get smaller and smaller in magnitude and eventually disappear in the limit. 

We are now in a position to formulate a rotation matrix whose elements 
correspond to the directional cosines of the x, y, and z axes of the rotated object. 

The rotation matrix = R^RpR^, 


where 


and 


Ra = 


cosa sina 0 
-sina cosa 0 , 
0 0 1 



1 0 0 
0 cosP sinP , 
0 -sinp cosP 


cosy 0 -siny 


Rt = 


0 1 
siny 0 


0 

cosy 
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Subsequently, 




coscxcosy - sinasinfisiny 
-cospsina 

sinycosa + cosysinasinP 


cosysina + sinysinPcosa 
cospcosa 

sinasiny-cosysinPcosa 


-sinycosP 

sinP 

cosPcosy 


where 


(3.40) 


a _ ^otj, P = ^Pj, and y = £y, . n corresponds to the iteration where all the rota- 
i=l i=l i=i 

tion terms go to zero in the limit. 

Once the rotation terms, i.e., xy, yz, and xz are eliminated, the three- 
dimensional surface has the representation of 


F(x,y,z) = Ax 2 + By 2 + Cz 2 + 2Px + 2Qy + 2Rz + D = 0, (3.41) 

where A, B, C, P, Q, and R are the coefficients resulting after the elimination of the 
rotation terms. A natural question to ask is: Can the terms of the first degree be elim- 
inated by means of a translation? The answer is sometimes they can and sometimes 
they cannot. The case, where the term can be eliminated, is supported by the follow- 
ing theorem. 

3.3.6 Translation of the Rotated Object 

Theorem 3.2. The terms of the first degree of an equation of a quadric surface 
can be eliminated by means of a translation if and only if the surface has a center, in 
which case the first degree terms are eliminated if and only if the new origin is a 
center [24] . 

The method of completing squares is the easiest to determine the coordinates of 
the new origin. Consider Equation (3.41). Grouping the like terms. 


=> 


Ax 2 +2Px +By 2 + 2Qy + Cz 2 + 2Rz + D = 0 
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r " 

x 2 + 2P-7- 
A 

+ B 

y 2 + 2Qf 

+ c 

^ + 2Ri 


+ D = 0 . 


Upon completing squares, we get 


r 

2 


2 

- ^ 

2 

P 

X + — 

+ B 

y + Q- 

+ C 

R 

z + — 

+ D - 

A 

W J 




C 

- 


P 2 Q 2 R 2 
A + B + C 


= 0 , ( 3 . 42 ) 


where -P/A, -Q/B, and -R/C are the coordinates of the new origin. 


3.4 Summary and Problem Identification 

All of the above procedures performed until now result in a second degree poly- 
nomial describing an unknown object, the center of the object lying at the origin of 
our coordinate system. Had the test data been simulated, the three-dimensional 
discriminant approach which was mentioned in Chapter One could be used to describe 
and recognize the object. Since the test data is not simulated, we should utilize a 
recognition algorithm which will distinguish and recognize each of the test surfaces 
from one another. 

The intersection of a surface with a plane generates a curve. The nature of this 
curve depends solely on what type of object is intersected and with which particular 
plane and in which orientation. Since we have no knowledge of the surface type, a 
priori, one approach is to intersect the surface with a series of planes. We need to 
determine the optimum number of planes which will uniquely characterize each of the 
quadric surfaces. 

Our goal is to derive a consistent method for determining the minimum number 
of planes necessary to intersect a given quadric surface so that the generated conics 
uniquely characterize the surface. This goal includes the derivation and formulation of 
the angular bounds for which a particular plane intersecting a surface generates the 
same two-dimensional curve. In summary, each of the quadric surfaces is represented 
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by a unique five-tuple, whose elements signify the presence or absence of the follow- 
ing curves: circle, ellipse, hyperbola, parabola, and a line. 

Chapter Four covers the description and recognition of each of the three- 
dimensional surfaces we have above mentioned in Section 3.2. A distinct pattern vec- 
tor is obtained for each of the surfaces. 



CHAPTER FOUR 


QUADRIC SURFACE CHARACTERIZATION AND RECOGNITION 
4.1 Introduction 

Our proposed method utilizes a two-dimensional discriminant which is a measure 
for distinguishing curves. Since the ten generated coefficients described in Section 
3.3.3 of Chapter Three give a three-dimensional representation of the surfaces, we pro- 
pose to identify the quadrics using the information resulting from the intersection of 
the surface with different planes. If the surface is one of those considered for the 
recognition process (see figures 3-1, 3-2, and 3-3), there are five possible two- 
dimensional curves that may result from such intersections: (i) a circle, (ii) an ellipse, 
(iii) a parabola, (iv) a hyperbola, and (v) a line. Thus, a feature or pattern vector with 
five independent components can be formed for characterizing each of the surfaces. 

The two-dimensional discriminant criteria we use to recognize each of the two- 
dimensional curves created by planes intersecting the various quadric surfaces is dis- 
cussed in Section 4.2. In Section 4.3 the results of Chapter Three are used to com- 
pletely implement our recognition algorithm. Concomitantly, we derive a consistent 
method for determining the minimum number of planes that are necessary to intersect 
a given three-dimensional surface so that the generated conics uniquely characterize 
the surface. The formulation of a three-dimensional discriminant similar to the two- 
dimensional discriminant is presented in Section 4.4. The mapping between the expli- 
cit and implicit representations of quadric surfaces is also examined in this section. 
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4.2 Two-Dimensional Discriminant 

Given a conic of the form 

F(x,y) = Ax 2 + Bxy + Cy 2 + Dx + Ey + F = 0, 

the discriminant 8 = B 2 - 4AC characterizes it as one of the following [30]: 

If 8 = B 2 - 4AC < 0, then the conic is an ellipse or a circle. 

If 8 = B 2 - 4AC = 0, then the conic is a parabola. 

If 8 = B 2 - 4AC > 0, then the conic is a hyperbola. 

Our objective is to derive a consistent method for determining the minimum 
number of planes required to intersect a given three-dimensional surface so that the 
generated conics uniquely characterize the surface. This includes the derivation and 
formulation of the angular bounds for which a particular intersecting plane yields the 
same two-dimensional curve. 

The three-dimensional surfaces (objects) to be recognized are listed below: 

(a) an ellipsoid, 

(b) a circular cylinder, 

(c) a sphere, 

(d) a quadric cone, 

(e) a hyperboloid of one sheet, 

(f) a hyperboloid of two sheets, 

(g) an elliptic paraboloid, 

(h) a hyperbolic cylinder, 

(i) a parabolic cylinder, 

(j) a hyperbolic paraboloid, and 

(k) a parallelepiped. 
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4.3 Quadric Surface Description and Representation 

As discussed in Section 3.2 of Chapter Three, we now assume that the three- 
dimensional objects have undergone two basic transformations, rotation and translation. 
Consequently the product terms in the representation F(x,y,z) for a particular surface 
have been eliminated and the center of the surface lies at the origin of our specified 
coordinate system. As illustrated in Figure 4-1, all of the surfaces are contained in the 
xy plane with their centers at O (the origin). For each surface, the characterization is 
performed in two steps. Initially we consider the intersection of each object with two 
planes (horizontal and vertical). This step does not require that the surface undergoes a 
translation transformation. We refer to plane 1 as the one that intersects the object 
parallel to the xy plane, i.e., z constant. Also refer to plane 2 as the one that inter- 
sects the object parallel to the xz plane, i.e., y constant. In the second step, the 
minimum set of intersecting planes needed to yield a unique feature vector (the various 
curves serve as features) is determined. In this step we assume that the object has 
undergone the translation transformation. The following sections describe the 
representation procedure for each of the quadric surfaces listed in Section 4.2. 

4.3.1 Ellipsoid 

Step 1: 

Consider the equation of an ellipsoid resting on a plane parallel to the xy plane 
and its axis of revolution parallel to the z axis. Equation (3.1) reduces to the form 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2px + 2qy + 2rz + d = 0, (4.1) 


which further reduces to 
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2 

■ 

2 


2 

x + £ 


y+£ 


r 

z + — 


a 


b 


c 


L J 



- + - 

_ 

— 


a 


- 1 = 0 , 


(4.2) 


2 2 j2 

where a>0, b>0, c>0 and we have assumed the scaling d = + 1 . 

a b c 

It should be noted that the coefficients a, b, c, p, q, r, and d are all known; 

and vy are the semi-major and minor axes of the ellipsoid, respectively; and 

[-p/a, -q/b, -r/c] are the coordinates of the center of the ellipsoid. 

Consider the intersection of the ellipsoid with plane 1, i.e., z = k, where 

— - '\Jl- < k < — + '\/^-, then, 

C 1 C c V c 


(y + y 


(x + £) 2 


a 


\ (ck + r) 2 

b 


be 


a 


(ck + Tf 
ac 


- 1 = 0 , 


(4.3) 


which is the equation of an ellipse. 

Let’s now consider the intersection of the ellipsoid with plane 2, i.e., y 

'\fr> then - 


= k. 


where -j3- ~ \j ~ 
b “ b 


<k< ^.+ 


(x + ^-) 2 (z + — ) 2 

a + C 1 = 0 , (4.4) 

1 (bk + q) 2 1 (bk + q) 2 

a ab c be 

which is again the equation of an ellipse. For the case when the two minor axes are 
equal, the surface is called a spheroid. Also, when all the axes are equal, i.e., a = b = 
c, the surface is a sphere. Intersection of the sphere with planes is discussed in Sec- 
tion 4.3.3. 
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Step 2: 

As mentioned before, we assume that the ellipsoid has undergone a translation, 
such that its center aligns with the origin of our desired coordinate system as shown in 
Figure 4-2. Hence its representation can be assumed as 



where A, B, and C are the major and minor semi-axes, respectively, of the ellipsoid. 
As seen in step 1, intersection of the ellipsoid with any Z=lkl,-C<k<C, will be 
an ellipse. Let us now determine the bounds within which inclined sub-planes of Z = 

I k I still result in an elliptic intersection with the ellipsoid. 

Consider the points E(A,0,0), F(0,B,0), and G(0,0,K), where K > 0. The equa- 
tion of the plane containing these points is: 

BKX + AKY + ABZ - ABK = 0. 

Solving for Z and substituting in Equation (4.5) yields the curve of intersection: 

X 2 (B 2 C 2 + B 2 K 2 ) + Y 2 (A 2 C 2 + A 2 K 2 ) + 2AK 2 BXY + ••• = (). (4.6) 

In the above equation only terms which are necessary to determine the intercepted 
curve are retained. Proceeding with the discriminant test, 

5 = 4A 2 B 2 [-C 4 - 2C 2 K 2 ]. 

Since the discriminant is always negative, the intercepts are ellipses. Angular bounds 
in terms of an angle are not needed in this case, since the only occasion the intercepts 
are different than ellipses is when two of the semi-axes are equal. Under that cir- 
cumstance, we arrive at a circular intercept. Figure 4-3 illustrates vertical planes inter- 
secting the ellipsoid. Table 4-1 summarizes the result obtained above. 
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DETAILED VIEW: HORIZONTAL INTERSECTIONS 


z 



Figure 4-2. The plane parallel to the x-axis and all its inclined sub-planes generate 
ellipses. In the case of a spheroid all intersections are ellipses except when the 
plane is parallel to one of the axes under which case the intersection is a circle. 
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DET AILED VIEW: VERTICAL INTERSECTIONS 


z 



Figure 4-3. The plane parallel to the z-axis and all its inclined sub-planes generate 
ellipses. In the case of a spheroid all intersections are ellipses except when the 
plane is parallel to one of the axes under which case the intersection is a circle. 
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PLANE 

INTERSECTION 

Z = K 

Ellipse 

Y = K 

Ellipse 

Any inclined sub-planes to Z=K, Y=K 

Ellipse 


Table 4-1. Intersection of ellipsoid with planes. 


4.3.2 Circular (elliptic) cylinder 
Step 1: 

Consider the general representation of a circular cylinder resting on a plane paral- 
lel to the xy plane and its axis of revolution parallel to the z axis. It’s representation 
then reduces to 


F(x,y,z) = bx 2 + by 2 + 2px + 2qy + d = 0, 


(4.7) 


which is the same as 


F(x,y,z) = 


- 

2 

■ 

2 

y + — 


x + £ 


3 b 


b 



— + - 

L J 

— 


_ 1 _ 

b 


b 


- 1 = 0 , 


2 2 

only if d = — — t- — 1 and also b > 0. 

b b 

In the case of the elliptic cylinder, Equation (4.7) becomes 
F(x,y,z) = ax 2 + by 2 + 2px + 2qy + d = 0, 


(4.8) 


which further reduces to 
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F(x,y,z) = 


* 

2 

“ " 

2 

y + -9- 


x + P 


7 b 


a 



- + -! 


— 


b 


1 = 0 , 


(4.9) 


2 2 

only if d = — 1 and a > 0, b > 0. 

b a 

Intersection of the circular or elliptic cylinder with plane 1 would not affect its 
representation, since it is independent of the variable z. Hence the resultant curve 
intercepted is the same as represented by Equations (4.7) or (4.9), which is an equation 
of a circle or an ellipse, respectively. 

Consider the case where the circular cylinder is intersected with plane 2, i.e., y — 
k, where - '\J^~ < k < Then, 


(4.10a) 


v-u P 

2 

1 

bk + q 

X H 

a 

" b ” 

b 

„ j 


Solving for x generates the equation of a pair of parallel lines. A similar result is 
obtained when the elliptic cylinder is intersected with plane 2, namely 


(4.10b) 


ol 

2 _ _1_ _ b^ 

bk + q 

x + — 
a 

a a 

b 


Step 2: 

As with the ellipsoid, consider the elliptic cylinder to have undergone the transla- 
tion transformation. Its center is aligned with the origin of the coordinate system as 
shown in Figure 4-4. Let the height of the cylinder be 2L. The representation of the 
elliptic cylinder can be assumed as 


xi + li = 1 

A 2 B 2 


(4.11) 
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where A and B are the major and minor semi-axes of the cylinder. Intersection of the 
cylinder with any plane Z = I k I, - L < k < L, is an ellipse. The angular bounds 
within which an inclined plane will still result in an elliptic intersection is determined 
next. 

Consider the intersection of the plane passing through the points E(A,0,0), 
F(-A,0,K), and G(0,-B,0) with the cylinder as shown in Figure 4-4. The equation of 
the plane containing these points is: 

BKX - AKY + 2ABZ - ABK = 0. 

Solving for X, 

„ AKY - 2ABZ + ABK 
X_ BK 

Substituting X in Equation (4.11) results in 

2K 2 Y 2 + 4B 2 Z 2 - 4BKYZ +... = 0. 

The discriminant results in a quantity less than zero. Hence the intersection is an 
ellipse. 

Given any two planes, ajx + b^y + CjZ + dj = 0, and a 2 X + b 2 y + ^2 Z + d 2 = 0, 
the angle of intersection is given as 

laja 2 + bjl >2 + ^1^2^ 

cos0 = , = ■ 

Va 2 + b, 2 + cj^/a 2 + b| + C 2 

Hence, in the above case the intersections with respect to the plane z = 0 and all the 
planes inclined to it (which we will refer to as inclined-sub planes), yield ellipses for 

2AB 

V(A 2 K 2 + 4A 2 B 2 + B 2 K 2 ) 


COS0 < 
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DETAILED VIEW: HORIZONTAL INTERSECTIONS 



Figure 4-4. Plane 1 and the planes parallel to it within the range -L to L 
Oength of the cylinder) intersect the cylinder in parallel lines. Plane 2 and 
plane 3 are the inclined sub-planes of plane 1 which determine the maxi- 
mum range or inclination (with plane 1) wherein similar curves ( ellipses) 
are generated . 0 is the angular bound for the inclination in terms of an 
angle. 
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The angular bounds with respect to the plane X = K -A < K < A, and its inclined 
sub-planes is determined next 

Intersection of the plane X = K or Y = K and the cylinder results in an inter- 
section of a pair of straight lines. The equation of the plane passing through the points 
H(0,B,-L), I(0,-B,-L), and J(K,0,L), I K I > 0, as shown in Figure 4-5, is 


KZ - 2LX + LK = 0. 


Solving for X, 


X = 


K(Z + L) 
2L 


Substituting in Equation (4.11), yields the interception 

K 2 (Z + L) 2 = j 

4L 2 A 2 B 2 


which is an ellipse. 

All intersections of the inclined plane X = K, I K I > 0, yield degenerate ellipses. 
In terms of the angle of intersection. 


cos0 < 


2L 

a/ A 2 + 4K 2 ’ 


Figure 4-6 illustrates a lateral view of all the possible curves intercepted by the inter- 
section of the cylinder and the planes. Table 4-2 summarizes the results obtained 


above. 
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DETAILED VIEW : VERTICAL INTERSECTIONS 



Figure 4-5. Plane 1 and the planes parallel to it within the range -a to a 
intersect the cylinder in parallel lines. Plane 2 and plane 3 are the inclined 
sub-planes of plane 1 which determine the maximum range or inclination 
(with plane 1) wherein similar curves (degenerate ellipses) are generated. 
0 is the angular bound for this inclination in terms of an angle. 



LATERAL VIEW 

INTERSECTION OF A PLANE AND A QUADRIC CYLINDER 



Figure 4-6. Plane PI and its inclined sub-plane generate ellipses. Though plane P2 
generates a pair of lines, its inclined sub-planes start generating degenerate ellipses 
as the inclination start to increase. 
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PLANE 

INTERSECTION 

Z = K 
X = K 

Inclined sub-planes of Z=K 
Inclined sub-planes of X=K 

Circle, Ellipse 
Lines 
Ellipse 
Lines 


Table 4-2. Intersection of quadric cylinder with planes. 


4.3.3 Sphere 
Step 1: 

As mentioned in Section 4.1., the sphere is a special case of an ellipsoid, 
the three semi axes are all equal. Equation (4.1) thus reduces to 

F(x,y,z) = ax 2 + ay 2 + az 2 + 2px + 2qy + 2rz + d = 0, 
which further reduces to 


a 


d 2 a 2 r 2 

only if d = — + — + 1. 

aba 



2 

- 

2 

' 

2 

X + -E- 


y + A 


r 

z + — 


a 


a 


a 



— + - 

_ 

— + - 

L J 

— 


a 


- 1 = 0 , 


Consider the case when the sphere is intersected with plane 1, i.e., z 

— - a/-^ < k < — + a/-^-. Then, 
a v a a ” a 


= k, 


where 


(4.12) 


(4.13) 


where 
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x + 


+ 5-1 

a 


L 

J 

T~ * 

2 A 


J_ _ 

ak + r 

i „ 

ak + r 

a 

a 

a 

a 

I*. * 


-1=0, 


(4.14) 


which is the equation of a circle. 

A similar equation results when the sphere is intersected with plane 2, in which 
case y = k, where — - *\j ~~ < ^ < + ^ an< ^ subsequently Equation 

(4.13) becomes 



X + 2 - 

a 

2 

i 


r 1 

ak + q 

2 

a 


a 





r 

z + — 
a 

2 

1 


r 1 

ak + q 

2 

a 


a 

J 



- 1 = 0 . 


(4.15) 


Step 2: 

Figure 4-7 illustrates the sphere which has undergone translation and has its 
center aligned with the origin of our desired coordinate system. The representation of 
the sphere thus becomes: 


A 2 


+ 



(4.16) 


where A is the radius of the sphere. As seen in step 1, intersection of the sphere with 
any Z = I K I, -A < K < A, will be a circle. Next, we determine the bounds within 
which inclined sub-planes of the Z = I K I plane still result in circular intersections 
with the sphere. 
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Consider the points E(0,0,K), F(A,0,0), and G(0,-A,0), where K > 0. The equa- 
tion of the plane passing through these points is 

-YK + AZ + KX - AK = 0. 

Solving for Z and substituting in Equation (4.16), yields the equation of the intercept 
as 

(A 2 + K 2 )X 2 + Y 2 (A 2 + K 2 ) - 2K 2 XY + • • • =0, 

where only the necessary terms to determine the nature of the intercepted curve are 
retained. Proceeding with the discriminant test, 

8 = -4A 2 [1 + 2K 2 ]. 

Since the discriminant is negative, the intercepts are ellipses or circles. Angular 
bounds are not needed since none of the other curves are ever intercepted. Similar 
results are obtained while considering inclined sub-planes of X = I K I or Y = I K I . 
Figure 4-7 shows a lateral view of all the curves intersected in a sphere by various 
planes. Table 4-3 summarizes the various results obtained above. 


PLANE 

INTERSECTION 

Z = K 

Circle 

Y = K 

Circle 

X =K 

Circle 

Any inclined sub-planes to X=K, Y=K, and Z=K 

Ellipse 


Table 4-3. Intersection of sphere with planes. 
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INTERSECTION OF A PLANE AND A SPHERE 



Figure 4-7. The intersection of plane and a sphere results in a circular 
line of intersection. 
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4.3.4 Quadric circular (elliptic) cone 
Step 1: 

The general representation of a circular cone on a plane parallel to the xy plane 
and its axis of revolution parallel to the z axis is 


F(x,y,z) = bx 2 + by 2 + cz 2 + 2px + 2qy + 2rz + d = 0, 

p 2 q 2 r 2 

where be < 0 and d = — + — h • 

b b c 

From Equation (4.17), upon completing squares, we have 


(4.17) 


F(x,y,z) = b 


- 

2 

■ 

2 

' 

x+£ 

+ b 

y + -J 

+ C 

r 

z + — 

b 


J b 


c 


+ 1=0. (4.18) 


2 2 J2 

Since d = — + — + — , Equation (4.18) becomes 
b b c 


F(x,y,z) = 


- 

2 

■ 

2 

' 

2 

x + ^ 




r 

Z H 


b 


} b 


C 


L .... J 

— 4- - 


— _ - 


— 


b 


b 


= 0 . 


In the case of the elliptic cone. Equation (4.17) reduces to 


(4.19) 


F(x,y,z) 


- 

2 

- 

2 

l 

2 

x + P 




z+ - 


a 


b 


c 



— + - 

L J 

! 

Z 

— 


a 


b 


= 0 , 


(4.20) 


where ab > 0, ac < 0, and be < 0. If c < 0, i.e., b > 0, the intersection of the cone 
represented by Equation (4.19) with plane 1, i.e., z = k, where — *\J ~ 

< k <— + 'V — — , would generate 
c V c 
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- 

2 

■ 

2 

' 

2 

x + £ 


y + r 


k + — 


b 


7 b 


c 




L J 

— = - 


t 


b b c 

where — — is a positive quantity. The above equation is that of a circle. 


(4.21) 
The elliptic 


cone on the other hand which is represented by Equation (4.20), upon intersection with 
plane 1, i.e., z = k, where < k ^ + V ’ WOuld generate 


r 

2 

- 

2 


2 

x + ^ 


y+£ 


k+ — 


a 


b 


c 


L - 

— + - 

_ 

— = - 

r ~ 



(4.22) 


a 


which is an ellipse. The intersection of the circular cone with plane 2, i.e., y-k, 

Vf 


where 

b 


< k < — + \J — , would generate 
b b V b 


r 

2 

- 

2 

' 

2 

x + ^ 


r 

z + — 




b 


c 


b J 



-1 1 


b c b 

where — — is a positive quantity. Equation (4.23) represents a hyperbola, 
c 

result is obtained when the elliptic cylinder is intersected with plane 2. 


(4.23) 
A similar 


Step 2: 

The quadric representation of the elliptic cone illustrated in Figure 4-8 is 

Y 2 Y 2 Z 2 

= 0 . 

A 2 B 2 C 2 


(4.24) 
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Intersection of the cone with horizontal planes z = k, where -c < k < c, generates 
ellipses as intercepts. Let us consider the horizontal plane Z = -C and determine the 
various intercepts formed by its inclined sub-planes. The equation of the plane passing 
through the points E(A,0,-C), F(0,-B,-C), and G(0,0,L) where -C < L < C, is 

-A(C+L)Y + ABZ + B(C+L)X - ABL = 0. 

Substituting Z in Equation (4.24) results in 

B 2 [C 2 - (C+L) 2 ]X 2 + A 2 [C 2 - (C+L) 2 ]Y 2 - 2AB(C+L) 2 XY + .... = 0, 

thereafter, 

5 = 4A 2 B 2 [(C+L) 4 - (L 2 +2LC) 2 ]. 

Analyzing 6 leads to the following bounds: 

For L > 0 the intersections are hyperbolas. 

C 

For all values of L, -C < L < O, except for L=-C+^=-, the intersections are 


ellipses. 

For the one particular case where L=— C+-^~ -, the intersection is a parabola. In 

terms of 0, the angle between the Z = -C plane and its inclined sub-plane is 

AB 

COS0 = , == = = = ~ - 

V(A 2 (C+L) 2 + A 2 B 2 f B 2 (C+L) 2 ) 

Next, consider the intersections formed by the plane X = 0 and its sub-planes. 
Substituting X = 0 in Equation (4.24) leads to the intersection 


Y^ 

B 2 



which is a degenerate hyperbola. For all -A < X < A, the intercepts are hyperbolas. 
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DETAILED VIEW : HORIZONTAL INTERSECTIONS 



Figure 4-8. Plane PI and planes parallel to it within the range -c to c (except 
the one passing through the origin) generate ellipses. Plane P2 is the inclined 
sub-plane which denotes the maximum inclination or range (of plane PI) wit- 
hin which ellipses are generated. 0 is the angular bound in terms of the angle. 
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The equation of the plane passing through the points H(0,-B,-C), I(0,B."C), and 
J(L,0,C), where L > 0 is 


LZ - 2CX + LC = 0. 


Solving for Z and substituting in Equation (4.24) leads to the representation of the 
intercept as 

X 2 [L 2 - 4A 2 ]B 2 + L 2 A 2 Y 2 + 4A 2 B 2 XL + • ■ * =0. (4.25) 


Solving L 2 -4A 2 indicates the following conditions for the various intercepts: 
For L = 2A, the intercept is a parabola. 

For all values of L, -2A < L < 2A, the intercepts are hyperbolas. 


For all L > 2A, the intercepts are ellipses. 

Figure 4-9 illustrates all of the above results. The angle between the X = 0 plane 
and its inclined sub-planes for the above obtained interceptions is 


cos0 < 


+2C 

V(L 2 + 4C 2 ) ' 


Figure 4-10 shows a lateral view of all possible curves intercepted in a quadric cone 
by the various planes. Table 4-4 summarizes all of the results obtained in this section. 


4.3.5 Hyperboloid of one sheet 


Step 1: 

The general representation of a hyperboloid of one sheet resting on a plane paral- 
lel to the xy plane and its axis of revolution parallel to the z axis is 

bx 2 + by 2 + cz 2 + 2px + 2qy + 2rz + d =0, (4.26) 

where base of the cylinder is circular and be < 0. 
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DETAILED VIEW : VERTICAL INTERSECTIONS 



Fieure 4-9. Plane PI and planes parallel to it within the range -b to b generates 
degenerate hyperbolas. Plane P2 is the inclined sub-plane which shows the outer 
region or the maximum inclination (of plane PI) within which hyperbolas are int- 
ercepted. 0is the angular bound in terms of the angle. Plane P3 is the only exce- 
ption where the intersection is a parabola. In this case the inclination of the plane 
P3 is equal to the base angle of the cone. 


LATERAL VIEW 


INTERSECTION OF A PLANE AND A CONE 



Figure 4-10. PI, P2, P3, and P4 are the four planes which generate all the 
intersections with the quadric cone. Plane PI which has the same base angle 
as that of the cone intercepts a parabola. Plane P2 intercepts a hyperbola. 
Plane P3 intercepts a circle and finally plane P4 intercepts an ellipse. (The 
quadric cone under question has a circular base). 



PLANE 

INTERSECTION 

Z = K 

Circle, Ellipse 

X = K 

Hyperbola 

Inclined sub-planes of Z=K, L>0 

Hyperbolas 

Inclined sub-planes of Z=K, -C < L < 0 

Ellipses 

c 

Inclined sub-planes of Z=K, L=-C+^=- 

Parabola 

Inclined sub-planes of X=K, L = 2A 

Parabola 

Inclined sub-planes of X=K, L < 2A 

Hyperbolas 

Inclined sub-planes of X=K, L > 2A 

Ellipses 


Table 4-4. Intersection of quadric cone with planes. 



74 


Equation (4.26) upon completion of squares reduces to 


F(x,y,z) = 


(x+ b )2 t (y+ b )2 , ~ <z + c )2 


b 


b 


-1 

c 


= 1, 


p 2 a 2 r 2 . 
where d = . + . — I A- 


If c < 0, i.e., b > 0, intersection of the hyperboloid with plane 1, i.e., z - k, where 
-r/c - *\J — - < k < -r/c results in 


L^2 


ty + fr 1 fr+ft 1 n (k+ t> 


b 


b 


-1 

c 


(4.27) 


where -1/c is a positive quantity. Equation (4.27) represents a circle. For a hyper- 
boloid with elliptic base, this intersection will be an ellipse. 

Intersection of the hyperboloid with plane 2, i.e., y = k, where -q/b - '\j — < 


k < -q/b +'\j generates 


(x + £) 2 (z + -) 2 (k + ^) 2 

b c | 


b 


-1 

c 


-1 

b 


where -1/c is a positive quantity. This equation is that of a hyperbola. Similar results 
are obtained when the hyperboloid has elliptic bases. 
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Step 2: 

As in the case of the other quadric surfaces, the elliptic hyperboloid of one sheet 
shown in Figure 4-11 is assumed to have undergone translation such that its center is 
aligned with the origin of the coordinate system. The axis of the hyperboloid coincides 
with the z axis. Under these conditions the quadric representation of the hyperboloid 

is 



The intersection of the hyperboloid with horizontal planes ranging from Z = 0 to Z 
= I K l are ellipses, where -C < K < C and A, B, and C are the semi-axes of the sur- 
face. The angular bounds of the various sub-planes with respect to the Z = 0 plane 
which intersects the hyperboloid in ellipses is determined next. 

As shown in Figure 4-11, the equation of the plane passing through the points 
D(A,0,0), E(0,-B,0), and F(K,0,C) where I K I > 0 is 

-ACY + B(A-K)Z - BCX - ABC = 0. 

Solving for Z and substituting in Equation (4.28) results in 

X 2 [B 2 (A-K) 2 - A 2 B 2 ] + Y 2 [A 2 (A-K) 2 - A 4 ] - 2A 3 BXY + • • • =0. 
Proceeding with the discriminant test, 

8 = 4A 2 B 2 [A 2 + K 2 - 2AK][A 2 - K 2 + 2AK], 

Since A and K are always positive, based upon the term 


[A 2 - K 2 + 2AK], 
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DETAILED VIEW : HORIZONTAL INTERSECTIONS 


z = o 



Figure 4-11. Plane 1 (z = 0) and all sub-planes parallel to it intersect the 
hyperboloid in ellipses. Plane 2 and plane 3 denote the maximum bound 
or inclination, within which the hyperboloid still intercepts ellipses. 
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a decision can be made whether the intersection is an ellipse, a hyperbola or a para- 
bola. Solving for K, we determine that for 

K = A(-V2 +1), the intersection is a parabola, 

K > A(-V2 + 1), the intersection is an ellipse, and 

K < A(-V2 + 1), the intersections are hyperbolas. The inclination of the plane at 
each of these intersections is given as 


B(A - K) 

/B 2 (A-K) 2 + P^C 1 + B 2 C 2 


Next, consider intersection of the plane Z = -C with the hyperboloid as shown 
in Figure 4-12. Substituting Z = -C in Equation (4.28) results in the intersection 



which is an ellipse as expected. To determine the bounds at which the inclined Z - 
-C plane still generates ellipses, consider the plane passing through the points 
G(L,0,-C), H(0,-B,-C), and I(M,0,K), where -C < K < C, I L I > I A L The equation of 

the plane results in 

B(C + K)X - L(C + K)Y - B(L - M)Z - (2BCL + BLK - BCM) = 0. 

Solving for Z and substituting in Equation (4.28) results in the intersection 

X 2 [C 2 B 2 (L - M) 2 - A 2 B 2 (C + K) 2 ] + Y 2 [C 2 B 2 (L - M) 2 - A 2 L 2 (C + K) 2 ] 

+ 2LBA 2 (C + K) 2 XY + • • • =0. 


Evaluating the discriminant leads to the following: 


5 = 4A 2 B 2 [L 2 A 2 (C + K) 4 - [C 2 (L - M) 2 - A 2 (C + K) 2 ][C 2 (L - M) 2 - L 2 (C + K) 2 ]] 
The bounds for the various intercepts are obtained as follows: 

M = L, the intersection is a parabola, 








M > L, the intersection is a hyperbola, and 
M < L, the intersection is an ellipse. 
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In terms of the angle, 


COS0 = 


B(L - M) 

Vb 2 (C + K) 2 + L 2 (C + K) 2 + B 2 (L - M) 2 


Next, consider the various intersections of the plane X - 0 and its inclined sub- 
planes as shown in Figure 4-13 with the hyperboloid. As seen before, for 
_ K < X < K, the intercepts are hyperbolas. The equation of the plane passing through 
the points J(0,B,-C), M(0,B,-C), and N(K,0,C) is 


-KZ + 2CX - KC = 0. 


Solving for Z and substituting in Equation (4.28) results in the intersection 

X 2 (K 2 - 4A 2 ) + A 2 K 2 Y 2 + • = 0- 

K. 


It is observed that for all K < I 2A I, the intersections are hyperbolas, 
case K = 2A, the intersection takes the form 

A 2 K 2 Y 2 + ^p~+ • • • =0, 

Iv 


However for the 


which is a parabola. Similarly for the case K > I 2A I the intersections are ellipses. In 
terms of the angle, the bounds for the plane X = 0 are 


COS0 < 


2C 

Vk 2 + 4C 2 


Figure 4-14 shows the lateral view of the various curves intercepted in a hyperboloid 
by various planes. Table 4-5 summarizes the results obtained in this section. 


80 


DETAILED VIEW : VERTICAL INTERSECTIONS 


z 4 


,\ plane 2 



Figure 4-13. Plane 1 and all sub-planes parallel to it intersect the hyperboloid in 
hyperbolas. The inclined sub-planes of plane 1 which are denoted in the above 
figure by plane 2 and plane 3 determine the maximum range or bound wherein 
hyperbolas are still intercepted. Beyond this range the hyperboloid intersects 
various planes in ellipses except the case when the plane makes an angle of 6, 
under which case the intercepted curve is a parabola. 
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LATERAL VIEW 

INTERSECTION OF A PLANE AND HYPERBOLOID OF ONE SHEET 






Figure 4-14. Plane PI intersects the hyperboloid in a parabola, plane P2 
and all planes parallel to it intersect the hyperboloid in hyperbolas. Plane 
P3 and all planes parallel to it in the range -c to +c intersect the hyperbo- 
loid in ellipses. 


C -Z. 
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PLANE 

INTERSECTION 

Z = K 

Circle, Ellipse 

X = K 

Hyperbola 

Inclined sub-planes of Z=0, K=A(-V2 + 1) 

Parabola 

Inclined sub-planes of Z=0, K>A(-V2 +1) 

Ellipse 

Inclined sub-planes of Z=0, K<A(-V2 +1) 

Hyperbola 

Inclined sub-planes of Z=-C, 1 Z 1 < C 

Ellipse 

Inclined sub-planes of Z=-C, 1 Z 1 > C 

Hyperbola 

Inclined sub-planes of X=K, K < 1 2A 1 

Hyperbola 

Inclined sub-planes of X=K, K = 2A 

Parabola 

Inclined sub-planes of X=K, K > 1 2A 1 

Ellipses 


Table 4-5. Intersection of hyperboloid of one sheet with planes. 


4.3.6 Hyperboloid of two sheets 
Step 1: 

Unlike the hyperboloid of one sheet, the hyperboloid of two sheets consists of 
two separate pieces. The quadric representation of a hyperboloid of two sheets lying 
on a plane parallel to the xy plane is 


bx 2 + by 2 + cz 2 + 2px + 2qy + 2rz + d = 0, (4.29) 

where the base of the hyperboloid is circular, be < 0. Completing squares results in 
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<*+$* (y + ij* 

F(x,y,z) = + j + 

"b b c 

where d = — h — + 1, and -1/c is a positive quantity. Intersection of the 

b b c 

object with the plane 1, i.e., z = k, where k I > 'J-Vc, results in 

(x + £) 2 (y + f) 2 (k+ 7 )2 

7 7 1+ zl ’ 

b b C 

where -1/c is a positive quantity. This equation is of a circle. For a hyperboloid with 
an elliptic base, this intersection will be an ellipse. However, when I k I = V(-l/cj, the 
intersection will result in a point. 

Consider the case when the object is intersected with the plane 2, i.e., y = k, 
where - q/b - Vi < k < -q/b + Vi- This intersection results in 


-<* + b )2 

_1_ 

b 


-(z + -^) 2 
c 

-1 

c 


< k+ b )2 

zL 

b 


+ 1 , 


which is an equation of a hyperbola. Similar results are obtained for a hyperboloid 
with an elliptic base. 


Step 2: 

As in the case of the hyperboloid of one sheet, the elliptic hyperboloid of two 
sheets is assumed to have undergone translation so that its center is aligned with the 
origin of the coordinate system as shown in Figure 4-15. The axis of the hyperboloid 
coincides with the z axis. Under these conditions the quadric representation of the 

hyperboloid is 
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(4.30) 


The intersection of the hyperboloid with the horizontal plane Z = K, 1 K I < C, is an 
imaginary ellipse. For I K I > C, the horizontal plane will intersect ellipses as seen 

from Equation (4.30). 

Consider the case where Z = -T, where T refers to the length of segment OG. 
Substituting in Equation (4.30) leads to the intersection 


X 2 Y 2 _ T 2 . 

a 2 t >2 /~<2 


which is an ellipse. Let us now determine the bounds wherein the inclined sub-planes 
of the plane Z = - T still intersect the hyperboloid in an ellipse. Equation of the 
plane passing through the points D(A,0,-T), E(0,-B,-T), and G(0,0,L), where -T < L < 
Tis 


A(T + L)Y - ABZ - B(T + L)X + ABL = 0. 
Solving for Z, and substituting in Equation (4.30) yields 


B 2 (C 2 - (T + L) 2 )X 2 + A 2 (C 2 - (T + L) 2 )Y 2 + 2AB(T + L) 2 XY+ • • • = 0. 
Discriminant 

8 = 4A 2 B 2 [(T+L) 4 - [(C 2 - (T+L) 2 ) 2 ]]. 

The bounds for the various curves are obtained as follows: 

For L = —=■ - T, the intersection is a parabola, 

<2 

p 

for L > -=■ - T, the intersection is a hyperbola, and 

yl2 


for L < -=• - T, the intersection is an ellipse. 

<2 
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DETAILED VIEW : HORIZONTAL INTERSECTIONS 
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Angle 0 at each of these values of L is determined as 

AB 

■ V-a 2 (-T4C) 2 + a 2 b 2 + B 2 (-T+C) 2 ' 

Next, consider the vertical plane X = 0 and its inclined sub-planes. Substituting X 
= 0 in Equation (4.30) yields the equation of a hyperbola. 

To determine the angular bounds of the various intercepts formed through the 
intersection of the inclined sub-planes and the hyperboloid, consider the plane shown 
in Figure 4-16. The equation of the plane passing through the points H(0,B,-T), 
I(0,-B,-T), and J(L,0,T) is 


LZ + 2TX + LT = 0. 


Solving for Z and substituting in Equation (4.30) yields 


X2 y2 _ [2TX + LT] 2 _ 


B 


2r>2 


L/C 


Expanding and re-arranging the terms, leads to the equation of the intercept as 

B 2 (L 2 C 2 - 4A 2 T 2 )X 2 + A 2 L 2 C 2 Y 2 + 4T 2 A 2 B 2 LX - A 2 B 2 L 2 (C 2 -T 2 ) + • • • =0. 

Based upon the term L 2 C 2 - 4A 2 T 2 which is the coefficient of X 2 , a decision can be 

made about the nature of the intercept. Since C < T, for all values L < 2A, the inter- 

2AT 

cept will be a hyperbola. The coefficient of X 2 will disappear when L = — , which 
case the intersection is a parabola. For all other cases, i.e., L > 2A, the intersections 
are ellipses. 

In terms of 0, where hyperbolas are intersected, the angle between the two planes 


is 


2T 

Vl 2 + 4T 2 


COS0 = 
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Figure 4-17 shows the lateral view of all the curves intercepted in a hyperboloid of 
two sheets by various planes. Table 4-6 sumrnarizes the results obtained above. 


PLANE 

INTERSECTION 

Z = K 

Circle, Ellipse 

Z = K, 1 K 1 > c 

Circle, ellipse 

Z = K, 1 K 1 = c 

Point 

Z = K, 1 K 1 < c 

Imaginary ellipse 

Z = -T 

Ellipse 

X = K 

Hyperbola 

Inclined sub-planes of Z=-T, L = -C 

Hyperbola 

c 

Inclined sub-planes of Z=-T, L = -^=- - T 

Parabola 

c 

Inclined sub-planes of Z=-T, L > -^=- - T 

Hyperbola 

c 

Inclined sub-planes of Z=-T, L < -^=- - T 

Ellipse 

Inclined sub-planes of X = K, L<2A 

Hyperbola 

2AT 

Inclined sub-planes of X - K, L< — 

Parabola 

Inclined sub-planes of X = K, L>2A 

Ellipse 


Table 4-6. Intersection of hyperboloid of two sheets with planes. 
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LATERAL VIEW 

INTERSECTION OF THE HYPERBOLOID OF TWO SHEETS 

WITH PLANES 



Figure 4 - 17 . Plane PI intersects the hyperbolid in a parabola, plane P3 and 
all planes parallel to it intersect the hyperboloid in ellipses, and finally, plane 
P2 and all planes parallel to it intersect the hyperboloid in hyperbolas. 
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4.3.7 Elliptic paraboloid 


Step 1: 

The quadric representation of the elliptic paraboloid resting on a plane parallel to 
the xy plane is 

ax 2 + by 2 + 2px + 2qy + 2rz + d = 0. (4.31) 

Equation (4.31) upon completing squares, reduces to the form 


(X + £) 2 (y + ■?->' z 

1_ + -t_ + 4-=0 


A, 2 


a 


b 


2r 


only if d = — I- . '\l ~t \l ~ are the semi-major and minor axes of the para- 

y a b ' V a ’ * b 

boloid, whereas l/2r is the height of the paraboloid. 

Consider the intersection of the elliptic paraboloid with the plane 1, i.e., z — k, 
where 0 < k < l/2r. The equation of the intercept is 


(x + ^) 2 
a 

l_ 

a 


+ ^ 2 
b 


-k 

_ 1 _ 

2r 


(4.32) 


where is a positive quantity. Equation (4.32) is that of an ellipse. 

l/2r 


-q 


Consider the intersection of the surface with the plane 2, i.e., y - k, where 
a / J_ < k < — — + \l — . The curve intercepted is the parabola 

y b b y b 


(x + ^) 2 
a 

1/a 


l/2r 


» + y 

l/b 
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Step 2: 

Unlike step 1, the elliptic paraboloid in this section has undergone translation. 
Hence, its center is aligned with the origin of the coordinate system as shown in Fig- 
ure 4-18. The axis of the paraboloid coincides with the z axis. Thus, the quadric 
representation of the surface is 


X 2 Y 2 

— + -i- + 2Z = 0. 
A 2 B 2 


(4.33) 


Intersection of the elliptic paraboloid with planes X = K and Y - K where 
< k < A, and — B <K < B, respectively, will yield parabolas as discussed in step 
1. Also, the planes Z = K, where K < 0, intersect the paraboloid in ellipses. Con- 
sider the intersection of the horizontal plane Z = -L (where L is the length of the seg- 
ment OG) and its inclined sub-planes with the paraboloid. The equation of the plane 
passing through the points D(A,0,-L) (where L is the length of the segment OG, and 
"a" is the semi-minor axis), E(0,-B,-L), and F(0,0,K), I K I > 0, is 

-A(K+L)Y + ABZ + B(L+K)X - ABK = 0. 

Solving for Z and substituting in Equation (4.32) yields the equation of the intercept 


as 


X 2 Yi _,A(K+L)Y-B(K+L)X + ABK 
A 2 + B 2 + AB 


Substituting K = -L will intercept the ellipse 


X 2 Y 2 01 
A 2 B 2 

Consider the case where K = 0. Under this condition the resultant intercept is 

ll + il-tf = 0 , 

A 2 B 2 B A 

which is an ellipse. Hence in the range -L < Z < 0, the inclined sub-planes of 
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DETAILED VIEW : HORIZONTAL INTERSECTIONS 



Figure 4-18. Plane PI and planes parallel to it intersect the paraboloid in 
ellipses. Plane P2 is one of the inclined sub-planes which determines the 
maximum inclination or range (of plane PI) within which ellipses are still 
generated. 0 is the angular bound in terms of the angle. 
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Z = -L, intercept ellipses. Analyzing the discriminant 8 leads to the following 
bounds: 

_ AB 

0 < COS0 < ^ a2(K+l) 2 _ A 2 B 2 _ B 2(L+K) 2 ) ' 

Next, consider the various intersections made by the plane X = 0 and its 
inclined sub-planes as illustrated in Figure 4-19. The plane X = 0 generates the 
intercept 

Y 2 

-±- + 2Z = 0, 

B 2 

which is a parabola. Consider the plane passing through the points H(0,-B,-L), 
I(0,B,-L), and J(N,0,M), where -L < M < 0, and 0 < N < A. The equation of the 

plane is found to be 

NZ - (L+M)X + LN = 0. 

Solving for X and substituting in Equation (4.32) yields the intercept 

N 2 Z 2 2LN 2 Z Yj _ 0 

(L+M) 2 A 2 A 2 (L+M) b 2 

which represents ellipses, except when N = 0. 

Hence all inclined sub-planes of the plane X = K, where -A < K < A, yield 
intercepts as ellipses. In terms of 0, 


cos0 < 


+(L+M) 

>/(N 2 - (L+M) 2 ) 


(4.34) 


Equation (4.34) denotes the angular bounds within which the intersections are all 
ellipses. Table 4-7 summarizes the various conics obtained when various planes inter- 
sect the elliptic paraboloid. 
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DETAILED VIEW : VERTICAL INTERSECTIONS 


t z 


P2 




Figure 4-19. Plane PI (x=0) and all planes parallel to it intersect the paraboloid in 
parabolas. The inclined plane P2 determines the range within which parabolas are 
still intersected. After an angular span of 0, the plane intersects the paraboloid in 
ellipses. 
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PLANE 

INTERSECTION 

X 

ll 

N 

Ellipse 

X 

ll 

X 

Parabola 

Y = K 

Parabola 

Inclined sub-planes of Z = - L 

Ellipse 

Inclined sub-planes of X = 0 

Ellipse 


Table 4-7. Intersection of elliptic paraboloid with planes. 


4.3.8 Hyperbolic paraboloid 
Step 1: 

Unlike the elliptic paraboloid, the hyperbolic paraboloid is symmetrical with 
respect to the xz plane, the yz plane and the z axis. Its representation is as follows: 

ax 2 + by 2 + 2px + 2qy + 2rz + d = 0. (4.35) 

In this case, however, ab < 0. Upon completing squares we have 

(x + -£-) 2 (y + -|j -) 2 

i b_ + _l_ = 0 

1 _j_ J_ 

a b 2r 

2 2 

only if d = and — is a positive quantity, 

abb 


Intersecting the surface with plane 1, 

i.e, z = k, results in 

(x + £) 2 
a 

+ ^ 2 - k 

j_ 

1 J_ ’ 

a 

b 2r 
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where -1/b is a positive quantity. This equation is that of a hyperbola. In the case 
when z = k = 0, it results in a pair of parallel lines which is a degenerate case of a 
hyperbola. 

Consider the case when the object is intersected with plane 2, i.e., y = k, then 

(x + — ) 2 (k + -3-) 2 

a , z a 

1 _L " _1 ’ 

a 2r b 

which is an equation of a parabola. The two planes considered in step I by them- 
selves prove sufficient enough to distinguish the hyperbolic paraboloid from all the 
other quadric surfaces considered for the recognition process. Hence, angular bounds 
to extract the regions where a unique set of features (curves) is obtained are not 
necessary in the case of this quadric surface. However, Figures 4-20 and 4-21 illus- 
trate the regions, if necessary, where extra features (curves) can be evaluated. Table 
4-8 summarizes the curves intercepted by planes 1 and 2 with the hyperbolic para- 
boloid. 


PLANE 

INTERSECTION 

Z = K 
X = K 
Y = K 

Hyperbola 

Parabola 

Parabola 


Table 4-8. Intersection of hyperbolic paraboloid with planes. 
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4.3.9 Hyperbolic cylinder 
Step 1: 

As in the case of a regular circular or elliptic cylinder, the surface of the hyper- 
bolic cylinder is parallel to the z axis. Subsequently, the variable z is not present in 
its quadric representation. It’s general representation when resting on a plane parallel 
to the xy plane is 

ax 2 + by 2 + 2px + 2qy + d = 0, (4.36) 

where ab < 0. Completing squares 

(X + £.)* (y + A) 2 

2 + 1=0 

j_ _! 

a b 

2 2 

only if d = — + — + 1. Also, -1/b is a positive quantity. Intersection of the 
a b 

cylinder with plane 1, i.e., z = k, generates a hyperbola. Since Equation (4.36) is 
independent of the variable z, the curve intercepted is the one represented by Equation 
(4.36). 

Intersection of the hyperbolic cylinder with plane 2, i.e., y = k results in the 
equation 


(x + -^-) 2 (k + ^-) 2 

a _ a _ . 

j_ ” 

a b 

which when solved results in a pair of straight lines. As in the case of the hyperbolic 
paraboloid, angular bounds to extract the regions where a unique set of features 
(curves) are determined are not necessary, since the two planes considered in step 1 
by themselves prove sufficient to distinguish this surface from all the other quadric 
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surfaces considered in the recognition process. A follow-up on the various inclined 

sub-planes of the z= k and y= k planes leads to a similar set of intercepts as with 

the x = k plane. Figures 4-22 and 4-23 illustrate the regions, if required, where extra 
features (curves) can be determined. Table 4-9 displays the intercepts formed when 
the hyperbolic cylinder is intersected with the two planes. 


PLANE 

INTERSECTION 

Z = K 
X = K 
Y = K 

Hyperbola 

Lines 

Lines 


Table 4-9. Intersection of the hyperbolic cylinder with planes. 


4.3.10 Parabolic cylinder 
Step 1: 

Unlike the two quadric cylinder considered before, i.e., the circular (elliptic) and 
the hyperbolic, this surface is parallel to the y axis. Hence the variable y is not 
present in its quadric representation. It’s general representation when resting on a 
plane parallel to the xy plane is 

f(x,y,z) = ax 2 + 2px + 2rz + d = 0. (4.37) 

Upon completing squares it reduces to 

(x + ^) 2 

a z 


1/a 


(4.38) 
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only if d = Intersection of the parabolic cylinder with plane 1, i.e., z - k, 0 < 

a 

k < 2r/ab, where b is any finite positive quantity signifying the width of the base of 
the cylinder, yields 

( * + a > 2 k 

I _L ’ 

a 2r 

which, when solved, results in a pair of parallel lines. 

Consider the intersection of the parabolic cylinder with plane 2, i.e., y = k. 
Since Equation (4.37) is independent of the variable y, the resultant curve intersected 
is the same as Equation (4.37), which is a equation of a parabola. As in the case of 
the hyperbolic paraboloid and the hyperbolic cylinder, angular bounds to extract the 
regions wherein a unique set of features (curves) are determined are is not necessary. 
The two planes considered in step 1 by themselves proved sufficient enough to distin- 
guish this surface from all the other quadric surfaces considered from the recognition 
process. Figures (4-24) and (4-25) illustrate the regions, if required, where extra 
features (curves) can be determined. Table 4-10 displays the intercepts formed when 
the parabolic cylinder is intersected with the two planes. 


PLANE 

INTERSECTION 

Z = K 
Y = K 

Lines 

Parabola 


Table 4-10. Intersection of the parabolic cylinder with planes. 
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DETAILED VIEW : VERTICAL INTERSECTIONS 



Figure 4-25. Planes parallel to the xz-plane (y=k) and all its inclined sub-planes 
intersect the parabolic cylinder in parabolas. 
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4.3.11 Parallelepiped 

Since planar surfaces cannot be represented with quadratic equations, we con- 
sider a plane of the parallelepiped. The general equation of a plane from Equation 
(4.1) is of the form 


2px + 2qy + 2rz + d = 0. (4.39) 

Intersection with plane 1, i.e., z = k, will generate 

2px + 2qy + d + 2rk = 0, 

which is the equation of a line. Similarly, intersection of the plane denoted by Equa- 
tion (4.39) with plane 2 will generate the line 

2px + 2rz + d + 2qk = 0. 

Table 4-11 summarizes the short results obtained for the parallelepiped. 


PLANE 

INTERSECTION 

Z = K 

Line 

X = K 

Line 

Y = K 

Line 


Table 4-11. Intersection of the parallelepiped with planes. 
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Table 4-12 summarizes the various curves (conics) derived from intersecting each 
of the eleven surfaces with the two planes z = k and y = k. These observations fol- 
low the results obtained in step 1 of each of the quadric surfaces. As seen from 
Table 4-12, the quadric cone and the hyperboloid of one and two sheets all generate 
similar curves. However, after using the results of step 2 (where angular bounds have 
been determined), we are able to distinguish each of the quadric surfaces from one 
another. Each of the quadric surfaces can be represented by a binary five-tuple, where 
the numeral 1 indicates the presence of a particular curve and the numeral 0 refers to 
the non-existence of that curve. Table 4-13 below presents the feature vector for each 
of the quadric surfaces. 

Quadric surfaces which seem to have identical feature vectors in the table above, 
get differentiated when the angular bounds theory as defined and derived for each of 
the surfaces (step 2) is applied. The next section briefly presents one other surface 
recognition approach which is very similar to the two-dimensional discriminant 
approach utilized to distinguish two-dimensional curves. It is one of our primary areas 
for future investigation. 

4.4 Mapping of Explicit to Implicit Representations for Quadric Surfaces 

Another objective which should be discussed is the formulation of a three- 
dimensional discriminant similar to the two-dimensional discriminant described earlier 
as means of recognizing three-dimensional objects. Consider the general quadratic 
representation of quadrics again, i.e., 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0. (4.40) 

This equation can be written implicitly, such that 


z = F(x,y) = Ax 2 + By 2 + Cxy + Dx + Ey + F. 


(4.41) 
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OBJECT 

PLANE 1 : x = k 

PLANE 2 : y = k 

Ellipsoid 

Ellipse 

Circle 

Circular cylinder 

Circle 

Line 

Sphere 

Circle 

Circle 

Quadric cone 

Circle 

Hyperbola, Parabola 

Hyperboloid of one sheet 

Circle 

Hyperbola, Parabola 

Hyperboloid of two sheets 

Circle, Point 

Hyperbola, Parabola 

Elliptic paraboloid 

Ellipse 

Parabola 

Hyperbolic cylinder 

Hyperbola 

Line 

Parabolic cylinder 

Line 

Parabola 

Hyperbolic paraboloid 

Hyperbola 

Line 

Parallelepiped 

Line 

Line 


Table 4-12. The various curves intercepted by the quadric surfaces when intersected 
with the planes z = k and y = k. 
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3-D SURFACE 

CIRCLE 

ELLIPSE 

PARABOLA 

HYPERBOLA 

LINE 

Ellipsoid 

1 

1 

0 

0 

0 

Circular cylinder 

1 

1 

0 

0 

1 

Sphere 

1 

0 

0 

0 

0 

Quadric cone 

1 

1 

1 

1 

1 

Hyperboloid of one sheet 

1 

1 

1 

1 

0 

Hyperboloid of two sheets 

1 

1 

1 

1 

0 

Elliptic paraboloid 

1 

1 

1 

0 

0 

Hyperbolic cylinder 

0 

0 

0 

1 

1 

Parabolic cylinder 

0 

0 

1 

o 

1 

Hyperbolic paraboloid 

0 

0 

1 

1 

1 

Parallelepiped 

0 

0 

0 

0 

1 


Table 4-13. Feature vectors (representing the prescence or absence of curves) for each 
of the quadric surfaces. 
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Equation (4-40) characterizes the complete surface through its representation, whereas 
Equation (4-41) characterizes surface patches on three-dimensional surfaces. As dis- 
cussed in Chapter Two, where surface curvatures are utilized to describe surface 
patches as being planar or curved (hyperbolic, parabolic, or elliptic), we wish to utilize 
the ten coefficients of Equation (4.40) in the form of a discriminant to represent 
patches on three-dimensional surfaces. In case we justify the existence of the implicit 
form, we would like to derive a mapping from F(x,y,z) to F(x,y); i.e., we would like to 
investigate the relations between A, B, C, D, E, and F and a, b, c, d, p, q, r and d. If 
this is possible, then we can attempt to derive a discriminant using A, B, C, D, E, and 
F, the combination of which can distinguish three-dimensional objects. 

We approach this problem in two directions. In the first approach we would 
numerically solve and derive relations for each coefficient in Equation (4-41) in terms 
of the coefficients of Equation (4-40). For example, while solving for F we arrive at its 
representation as 

^ -2r ± V4I- 2 - 4cd 

F= 2c ' 

Similarly, expressions for the coefficient B have been found to be 

„ . V2 (f+g+r) 2 - 4c(a+b+2h+2p+2q+d) 

B “ ± 4c ' 

Each of the above coefficients were derived while setting the variables x and y as 
zero. Similarly the remaining coefficients can be derived after solving several linear 
and non-linear equations. In the second approach, derivatives are utilized to obtain a 
pattern vector based on the coefficients of Equation (4-40). Rewriting Equation (4-40) 
in terms of a quadratic of z, 

-[cz 2 + 2fyz + 2rz + 2gxz] = f(x,y) = ax 2 + by 2 + 2hxy + 2px + 2qy + d, (4.42) 


Ill 


which is similar to Equation (4-41), i.e., to 


z = F(x,y) = Ax 2 + By 2 + Cxy + Dx + Ey + F. 

Differentiating Equation (4-42) with respect to each of the variables, x, y, and z, 
yields the following equations: 


— = 2fz = 2by + 2hx + 2q, 
dy 


and 


— = 2gz = 2ax + 2hy + 2p, 
dx 


— — = 2cz + 2fy + 2r + 2gx = 0. 
dz 

Each of these expressions are utilized individually in Equation (4-40) to yield an 
expression of the form 


Ax 2 + C 2 + Bxy + Ex + Fy + D = 0, 

from which the discriminant B 2 - 4AC again produces results which are either less 
than zero, equal to zero, or greater than zero. A list of pattern vectors which seem to 
be invariant has been derived for some of the quadric surfaces. Much more work has 
to be done on simulated as well as real data before arriving at definite conclusions. 

The theory developed in chapters Three and Four were experimented with several 
sets of real and simulated range data. Results of which are summarized in Chapter 


Five. 


CHAPTER FIVE 


EXPERIMENTAL RESULTS 


5.1 Introduction 

Section 5.2 discusses our study of median filtering on range images. Section 5.3 
explains the process whereby filtered range images of spheres, cylinders, and quadric 
cones undergo the recognition criterion. Subsequently, Section 5.4 discusses the appli- 
cation of the rotation alignment algorithm to the processed as well as simulated range 
images. Section 5.5 briefly presents the results obtained while using the three- 
dimensional discriminant approach to simulated and real range data. 

Range data obtained using a laser radar three-dimensional vision system is similar 
to intensity images obtained from a normal camera. However, instead of intensity 
(brightness) information, range (depth) information is available. Thus it is possible to 
interpret intensity information as range information when a range image is displayed 
on a image processing monitor. The nearer the object, the brighter it appears on the 
screen. 

The experimental work was performed in the following order : 

(i) The effect of median filtering on range images was studied. 

(ii) The proposed recognition scheme was applied to filtered range images. 

(iii) The quadric alignment algorithm was applied to simulated and real data. 

(iv) The three-dimensional discriminant approach was tested with simulated data. 
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5.2 Median Filtering on Range Images 

Range images of objects like spheres, cylinders and cones were segmented in 
order to separate the object from its background. The resulting image, which is 
referred to as the raw image, was then median filtered with mask sizes 3 x 3, 5 x 5 
and 7x7. 

Consider Figure 5-1 which is the actual range image of a sphere with its back- 
ground. Figure 5-2 is the image after segmentation. The effect of median filtering on 
Figure 5-2 can be observed in Figure 5-3 (3x3 mask). Figure 5-4 (5x5 mask), and 
Figure 5-5 (7x7 mask). The curvature sign map, which is discussed in Chapter 
Three, was then utilized to study the effect of median filtering on the original range 
image shown in Figure 5-2. Evaluating the first and second derivatives with respect to 
the x and y axes and comparing each of these maps determines whether or not the 
median filtering has altered the original range data to any extent. Figures 5-6a, 5-6b, 
5-6c, and 5-6d are the first and second derivatives with respect to the x and y axes, 
respectively, for figure 5-2. Similarly figures 5-7a, 5-7b, 5-7c, 5-7d, figures 5-8a, 5- 
8b, 5-8c, 5-8d; and figures 5-9a, 5-9b, 5-9c, 5-9d are the first and second derivatives 
for the images in figures 5-3, 5-4, and 5-5, respectively. 

In all of these figures, "+" is assigned to a particular pixel position if the magni- 
tude of the derivative (first or second) of that pixel is greater than the magnitude of the 
derivative (first or second) of the pixel to its right. Similarly is assigned to a par- 
ticular pixel position if the magnitude of the derivative (first or second) of that pixel is 
less than the magnitude of the derivative (first or second) of the pixel to its right. In 
the case when the magnitudes of the derivatives (first or second) of either pixels is the 
same, a blank is assigned. 

Sign maps are also generated to check the integrity of the image data before and 
after the filtering process. Based on the magnitude of the depth value of a pixel and 
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Figure 5-1. Raw range image of the sphere with its background. 



Figure 5-2. Range image of the sphere after segmentation. 
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Figure 5-3. 3x3 median filtered image of the raw sphere. 



Figure 5-4. 5x5 median filtered image of the raw sphere. 



Figure 5-5. 7x7 median filtered image of the raw sphere 
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Figure 5-6(a). First derivative w.r.t x-axis of the original sphere. 
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Figure 5-6b. First derivative w.r.t y-axis of the original sphere. 
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Figure 5-6c. Second derivative w.r.t x-axis of the original sphere. 
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Figure 5-6d. Second derivative w.r.t y-axis of the original sphere. 
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Figure 5-7a. First derivative w.r.t x-axis of the sphere filtered with a mask size 
of 3 X 3. 
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Figure 5-7b. First derivative w.r.t y-axis of the sphere filtered with a 
mask size of 3 X 3. 
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Figure 5-7c. Second derivative w.r.t x-axis of the sphere filtered with a 
mask size of 3 X 3. 
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Figure 5-7d. Second derivative w.r.t y-axis of the sphere filtered with a 
mask size of 3 X 3. 


125 


-+ 
+ + 
+ + 
-+ 
-+ 


- + + +-+-++- 

-++++-+4+-+++ 
- + + +4-4+4 


4 

4 


-+-+++-+-+-++++++“+ + 4 + 44 --++ 

4- + + + +++++++ +“++++++“++ 

+++ 4-4- 4-4- 

4- 4-4-4- 4- - 4-4- “4 4-4- 


-++-- 

+-+-- 

4 + 

+ 

- + 

+-+-“ 

4 + 4 

+ 

-4 

+ 

+4 + 

- 44-44 

. + 

+ 

--+-++ 

- 44 + 4 + 

— 

+ + 

+ 

44 + 4 + 4 


4-4- 4- -4- 4- 4- 4- 4- 4- 4- 


4 - 4 - 4 - 4 - - 4 - - 4 - 4 - 4 - 4 - 4 - 4 - 4 - 4 - - 4 - 4 - 4 - - 4 - 4 - 4 - - 4 - 4 - “ 4 - 4 - 4 - 4 - 4 - 4 - 4 - 4 - 4 - 4 4 - 4 - 4 - 4 - 4 - 4 - 

+_+++++ + 4-4 + 4+44 4 --4 4 -“4 -4 44 “444 + 4 4 4 444 - + + + + + + + + + + - + -4 + 44 + 4 + + 

+ - + + + + + + + ++ 44 - 4 -+- + + “ - 4 -4 + 4-4 + 4 + + + + - 4 - + +- + + + + +- + ++ 4 - 4 -- + 44 + + 44 

4 -+++++++ 4+44 +++ ++++++-++++-+++++++++++++++++ + ++ + “ ++++ + 4 + 

._+_+_++ 4444 + + — + 4 — + +4 + +— + + + + + + + 4 - + + + 4-4 + + + + + 4 + 4-4 + 4 + 


+-4-4+4+44- 

+ + + + + + + +++--+ + + + + + + + + + + + + + 4- + + + + 44 + -4-444 + + + + 44 + 4 

+ +--4-4 + 4 + + + 

4 + -- + +- 

+ 4-++-4-+ +--+ 4-4 + + + + + + + + + + + + + + + +- 

+4+4 -4+ -+++ ++ +--+++ +-++++-++++++++++++++++++++++++4+4+++" 

+--+-++ +- + ++ + + + + + +- + +- - + -- 44-4- + + - 

+ 4-44 4 + + + + + - + +- +4-4 + 4-+ + + + - 

+ +-++++ + + — + F +— 4 F + + ++— + + +— + + + + + + 4 + ++ 4' + + 4- — 

+ + --+ +--+ 4-- 

+--+--+--+ 

+- + + +--+--+ 

+ -4 + +4 + + + + + + + +-- + + + 4- - + 

4+-+-++444++++4++4+4++-++++4++++++++++++++++++ - +++++~+++++444++ 

+ 4 — + — + + 4 1-+++++ + + + - + + +“ + + + + + 4 4 + + + + + + + + + + + + + + + + + 4' + + + + *’ + + + 44 I-+ + + 

- + + + -4 + 4-4- + + + - + --+— + + + + 4 + 4 + -4- + + + -4 +-+++-+-- + + + - 

+-_+-_++-+- - + + + + - _+ + +4 + + + +- + +4 + + + + +-++ + - 

++-++--+--+ +--+--++-+--++++ +++++++ -+ ++ 4+ 4+ +4 4+ 4+ 4- 

4+ - + + + + --+ 4-+-- + --+--+ + 4- 4+ 4 + 44 + 44-4 + 4+ 4- 

+ ++ + + + + + + 4 + +- + + 4 + 4- ++- + 

+ + + - + --+-- + + + + + 4 -4 + ++ + + + 4 + 4+ 4+ “+-4- + 4 4 4- + -- + 

+ + -++-+-- + + + + + + + 4+ h++ + 4 + + 4 + + + +4 + 4 + 4-+— +— 4 + +“4 + 4 — 4 4 

+ + + + 4 4 4 44-4-4 + + + + + - 

+ + 4 + +““ 

++++-+444 + -++ +444++++ + ++++++ + ++++++ ++ ++++ 4+44++++ 4+ ++ 4+44+4 ++ ++ 
+++++44+4+++++4+-4++++++++++4+++++++++++++++++++++++++++++ +++++ 

+ +-+ 4 4 + + +-4 

+ + -+ - + + + + 4 +-- + --+ + “4 

+++++++++ +4+4+4 4+ --+-++ --+++- -++ -44 -4 -4+ -44+4 + 4+ -44+-++ +4+ 

+++4++-++4+4+4++++++++4--+4+++— ++++— 4+4— +—++++“-++++ 4+4+444+ 

+ + + 44++ + + 4444+ + + +4+ + + 4 + + 4 + + +4 + — + + + + + — + + h - + ++— 4 >- + + + + + -- +- 4 

+ + + 44 + + + ++ + 44 + + + + + 4 + + 44 + 4 + + + + + + 4 + + 4+ h + 4 + + — + + + + +4~+ + + -I +- + 


4 


4 


Figure 5-8a. First derivative w.r.t x-axis of the sphere filtered with a 
mask size of 5 X 5. 
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Figure 5-8b. First derivative w.r.t y-axis of the sphere filtered with a 
mask size of 5 x 5. 
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Figure 5-8c. Second derivative w.r.t x-axis of the sphere filtered with a 
mask size of 5 X 5. 


128 


++++++++++ 
--++++++++++ 
++++--++ + 


-+++++++++++++++++ ++ + 

++++++-++-+++++--+++++-+--+“++ ++++++ 

+ + — f- + + ++++~++ 

--+ +-+--+ + +- + -“+ ++ + 

++-4+-+-+ ++ ++ ++ +++++++- + 

__ + __ + ++ ++ +++ + +++--+++++++++-+++++-+++++++++++ ++++++++++ 

+++ +++ + + + + + - ++++++++++++ - --+++- - ” +" + 

h + ++“++“++++ 

- + ++--+--+--+-++++++ 

- + + + - “ - + + 

+ + + _ 4 + - + -i 1- + + + +++++++ 

-+++++++ + + ++ + + + + + " + ++++ “++++++++++ 


++++++++ 

- + + 


+ +- + 
+ 


+ + 
+ + 



++++++++ ' 




++++++++--+++++ + 

+ + --+ + +“+- 

+ ++++++ +-++++++““++“+" 

-++++++++++++ -+++++ ++++++++++++++++++++' 

4-4-44 + + + + + ++++++ + + + + + + + + + + + + - 

+ + + + + + + ++ + + + + + + + + + + + + + + + + + + ■! 

+4+444+ 

+ + + ++++++++++++H 

+ 4 * + + +”"' — 

4 44444 + 44444 + 4 + + + + + + + + - + +++ + + + + + + + + + + + + + + + + + + + + + + + + '*''*"*' + + + 4 ’' f ’ 

44444+ 4-4-4- + + + + + + + + + + + + + + + + + + + + + 

-+++++++-+-++++++++++++ ++ ++++++ -++ -++-+ 

-+4-4- + + ++ + — + + + + + + + + + + + + + + + + + + ++ + + + + + 

4.4-44.4-4- + + + + + + + + + + “ + + + ” 

4 444444 44444444+ 

h + 44444 4 + + + + + + + + ++ + + + + + + + + + 1-444444 + + + + + + + + + + ++ + + + 

4444444++++44++++++++++++++++++++++++++++++ ++++++++++++++++++ 

+ + + + + + 


++++++++++-+++-++ 


+ + - 


^++ + - 






-++ -+++-- 


4 444-4444444-4444444 + + - + + + + + + + +' 



4+44444+4444+++++++++++++ +++++++++++++++++++++++++++++++++ 

44444+-+ 444 + 44444 44+4 + 4 ++ + + - + 

-4 + + 4444--+- + + 4 + 44+ + + + + + + +++ + 

44+444 4+ ++++++++4+ +++++++ +4+4+++++++++++ 

+ + + + + + + + + + + + + + + + + + + 

+ -4 + + + + ++ - 

+ + + 44 “ “ “ “ ” 

++++++++ ++++++++++++++++++++++++++++++++++++ ++++++++++++++ 

+++ __ + ++++++ + +++++++++ +++++++++++ 4 -++++++++++++++++++++ +++ + 

44 + + 

4444 + 4-4 + 4 + 44 + 4 4 + + 4 + + - + 4 4 4 4 + + - + + ++--+++++++++“" 

4444+4444+44444+4-++++ +++444+++--++++++++++++++++++++++ 

+ 4- + + 4 + 4 + + + + + 4 + + + + + + 

4 + 4444 + + + + + ++ + + + + + + + + + + 4+ + — + +- + + + + + + + + + + + + + + + + + 

-++44 444 + 4+ 444 + 4 + + ++ + +--+ + +- + '■ 

+++++++++++++++++++++ 

444 


44 + -- + + + + 


- + + + 


Figure 5-8d. Second derivative w.r.t y-axis of the sphere filtered with a 
mask size of 5 X 5. 
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Figure 5-9a. First derivative w.r.t x-axis of a sphere filtered with a 
mask size of 7 x 7. 
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Figure 5-9b. First derivative w.r.t y-axis of a sphere filtered with a 
mask size of 7 x 7. 
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Figure 5-9c. Second derivative w.r.t x-axis of the sphere filtered with a 
mask size of 7 x 7. 


132 


-++++H 


-+ + + H 


- + 


H++++++++++++ 
- + + + + + -+ 


++++-+++-+“+-+++--■ 

+ -+“ + 

+ + 


++++-++++++++++“++++++++++ 

+ + 

+ -++ “+ + + + - + 

-+ + + + “+ + 

-++++- “ + - 

+ + + + + + + + + 

+ + +- + + +- 4 *+ + + + +~ + + + 

-+++++ -- + + + + +“" " + + + “ +_+ 

+ + + + + “+"++“ + ~ + 

+ + + - - -+ - ++ + ++ + 

+ + + + + + + + + - + + - 

- + _ - + + + - + +4.4-4- + + + + + -+- + + + 

+++-++--+ + + + + +++++ + 


+ + + 
+ - 


-- ~ + 
+ + + 

+ + 


+ 

+ 


++ + ++++ +++ + + 

__ _ + + + 

+ + + + 

+ ++ +-+ “+ 

------ ++ -+ ++ 

+ + +— +— +— + + + “+ — — — 

+-++ ++++++ ++++++++ +--++ -+++++-++-++++++--+++++++++++++++“ 

+++++++ ++++++++++++++++++++++++++++++++++++++++ ++++++++++++ 

+ + + + + + + + + + + ’ 

++++++++ ++++++++++++++++++++++++++' 


-++-++++++++++++++++ 


-+++++++++++++++++++++++++++++++++++ + 
++++-++++++++++++-++- 
++++-++++--+--++ - 


+ + + ++ 

++++++ 

+ + + + + + + + + + + + + + + 

+ + +— + + + + + + + + + + + + + 

+ + + + ++ - + + + + ++ + + + ++ + + + + + 

(. + + + -- + ++ + + ++ + + + + + + +- 

h++ ++++++++++ ++++-+ +++++++++++--+++++++++ 

+++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++ 

--+++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++ 

-++++++ 

+ + + + 


-++ - 


+ + - 


+ + + + + + + + + + + +++-++++-+-++--+ +- + - + + +-"- 

+ + 

-++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++ 

+ + + + + + + + ++ + + ++ + + + + ++ + 

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 

+++++++++++++++++++++++++ +++++++++++++++ 

+ + + + + + + + + + + + + + + + + + + 

+ + + + + + + + + + + + + 

+ + + + + 

++ ++ ++++++++++++++++++++++++++++++++++ ++++++++++++++++++++ 

+ + + + + + + + + + + + + + ++ + + + + + ++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 

++ — — — — — + + 

+ + + + + + + + + + +- + + + + + + + + + + + + + +-++ + + + + ++ + +- + ++ + + +- + + + + + + + + +“- 

+ + + + ++++++ ++++-++++ +++++++++--++++++++++++++++++++++ 

+ + +-+++ + + + + + + + + + + + + + + 

+ + ++ + + + + + +- + + + + + + + + + + +- - + + + -+ +-+ + + +++ + + + +++ + + + + “ 

+++++--+-++++++ + + + + + + + ++++“- + + + ' + + + ++ + + + + + + + 

+ + + + + + + + + + +++ + + + + + + + + + + + + + + + + + + + + + + + + +“+ + + +++ + + + + 

+ + + + -- + +++-- 


- + + + 


Figure 5-9d. Second derivative w.r.t y-axis of the sphere filtered with a 
mask size of 7 x 7. 
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its adjacent neighbor, a "+" or or " " (blank) is assigned to the pixel location in the 
sign map. Figure 5-10 is the sign map generated for the original raw image data of the 
sphere. Similarly figures 5-11, 5-12, and 5-13 are the sign maps generated for the 3 x 
3, the 5x5, and the 7x7 filtered images of the sphere. A careful observation of 
these sign maps suggests that only a small variation has been brought about due to the 
filtering process. 

The prime objective of median filtering is to remove salt and pepper noise in the 
range images and thus present a noise free range image for the evaluation of the 
objects coefficients [27]. It can be seen from figures 5-3, 5-4, and 5-5 that these filters 
met the objective. However, looking at the curvature maps it is observed that as the 
filter size increases, the apparent curvature is distorted relative to the original curva- 
ture. The 3x3 filtered image, being the closest to the original raw image, can be 
utilized for further processing and for describing the surface features. The validity of 
the curvature map calculations were checked using a "best fit analysis. 

Once the data files were obtained for each of the filtered images, the depth infor- 
mation of each of these files was converted into rectangular coordinates. The opera- 
tion manual for the laser radar three-dimensional vision system [31] describes the 
equations used for the transformations of the range information from spherical coordi- 
nates to rectangular coordinates: 


X = (R - L)sin0 f , 


(5.1) 


Y = (R - 


_ 5 _ 

cos 


0|— L)sin0 g cos0f, 


(5.2) 


and 


Z = (R - 


cos 


0(— L)cos0 g cos0f, 


(5.3) 


where 0 f is the horizontal scanning angle and 0 g is the vertical scanning angle. 
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Figure 5-10. Sign map generated for the original raw image of the sphere taking into 
consideration the magnitude of the depth value at a particular pixel and its neighboring 

pixel. 
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Figure 5-11. Sign map generated for the 3x3 filtered image of the sphere taking 
into consideration the magnitude of the depth value at a particular pixel and its 
neighboring pixel. 


Figure 5-12. Sign map generated for the 5 x 5 filtered image of the sphere taking 
into consideration the magnitude of the depth value at a particular pixel and its 
neighboring pixel. 


Figure 5-13. Sign map generated for the 7x7 filtered image of the sphere taking 
into consideration the magnitude of the depth value at a particular pixel and its 
neighboring pixel. 
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0 f = 25° - (horizontal pixel #) (0.1961 deg/pixel). (5.4) 

0 g = (vertical pixel #) (0.1961 deg/pixel) -25°. (5.5) 

L = 0.362m. (5-6) 

R is Range in meters = (0.00459 m/pixel)(Range pixel) + (n — 1/2), (5.7) 


where n is the electronic range in meters set by the operator. The cartesian coordi- 
nate information was then utilized for determining the coefficients which describe each 
of the three-dimensional surfaces. 

Experiments were conducted on range data for spheres and cylinders. Results of 
median filtering for one such set of range data is presented. 

Figure 5-14 is the actual range image of a cylinder. Figure 5-15 is the range 
image after segmentation. Similarly, figures 5-16 and 5-17 illustrate the 3x3 and the 
5x5 median filtered images of the cylinder range data. 

Curvature maps for studying the effect of median filtering on the range data are illus- 
trated in figures 5-18(a,b,c,d), figures 5-19(a,b,c,d), and figures 5-20(a,b,c,d), which are 
the first and second derivatives with respect to the x and y axes for the original 
cylinder image, the 3x3 filtered cylinder image, and the 5x5 filtered cylinder 
image, respectively. 

Sign maps similar to the ones derived for the sphere are generated for the 
cylinder and are shown in figures (5-21), (5-22), and (5-23). The figures correspond to 
the original, the 3x3 filtered image, and the 5x5 filtered images of the cylinder, 
respectively. Analyzing the curvature maps for the cylinder indicates the filtering pro- 
cess removed the noise and smoothed the image data without effecting significant dis- 
tortions. The sign maps, much like the curvature maps, seem not much affected by the 
filtering process other than some information (range data) being lost at the edges. 
Listed in tables 5-1 and 5-2 are the coefficients obtained for the original range images. 
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Figure 5-14. Raw range image of the cylinder with its background. 



Figure 5-15. Range image of the cylinder after segmentation. 
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Figure 5-16. 3x3 median filtered image of the raw cylinder. 



Figure 5-17. 5x5 median filtered image of the raw cylinder. 
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Figure 5-18a. First derivative w.r.t x-axis of the original cylinder. 
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Figure 5-18b. First derivative w.r.t y-axis of the original cylinder. 
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Figure 5-18c. Second derivative w.r.t x-axis of the original cylinder. 
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Figure 5-18d. Second derivative w.r.t y-axis of the original cylinder. 
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Figure 5- 19a. First derivative w.r.t x-axis of the cylinder filtered with a 
mask size of 3 X 3. 
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Figure 5- 19b. First derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 3 X 3. 
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Figure 5- 19c. Second derivative w.r.t x-axis of the cylinder filtered with a 
mask size of 3 X 3. 
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Figure 5-19d. Second derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 3 X 3. 
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Figure 5-20a. First derivative w.r.t x-axis of the cylinder filtered with a 
mask size of 5 X 5. 
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Figure 5-20b. First derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 5 X 5. 
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Figure 5-20c. Second derivative w.r.t x-axis of the cylinder filtered 
with a mask size 5X5. 
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Figure 5-20d. Second derivative w.r.t y-axis of the cylinder filtered with 
a mask size of 5 X 5. 
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Figure 5-21. Sign plot for the original cylinder. The sign M +" or is 
assigned depending whether the adjacent pixel has a range value lesser 
or greater than the pixel to its left. 
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Figure 5-22. Sign plot for the 3 x 3 filtered image of the cylinder. The sign 
"+" or is assigned depending whether the adjacent pixel has a range value 
lesser or greater than the pixel to its left. 
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Figure 5-23. Sign plot for the 5 x 5 filtered image of the cylinder. The sign 
"+" or is assigned depending whether the adjacent pixel has a range value 
lesser or greater than the pixel to its left. 
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the 3x3 filtered images, the 5x5 filtered images and finally the 7 x 7 filtered 
images of a sphere and a cylinder. 

These tables show that none of the coefficient sets describe a real sphere or 
cylinder with any certainty. The following procedure was utilized to determine which 
particular set of coefficients best describes the original range data of the object. A 
small surface patch of the object is chosen. In the quadratic form, 

F(x,y,z) = ax 2 + by 2 + cz 2 + 2fyz + 2gzx + 2hxy + 2px + 2qy + 2rz + d = 0, 

the coefficients a, b, c, d, f, g, h, p, q, and r are inserted and for each (x,y,z) of the 
object patch, the error is evaluated for each set of coefficients. A plot is generated in 
which every point of the surface patch is replaced with the numerals 1, 3, 5, and 7 
signifying that the minimum error was obtained for that particular set of coefficients. 
Figure 5-24 is one such plot for the raw segmented image of the sphere. Numeral 1 
refers to the situation when the original set of coefficients fits best, and similarly 
numerals 3, 5, and 7 are used depending whether the 3x3 or the 5x5 or the 7 x 
7 set of coefficients, listed in Table 5-1, of the sphere fits best. As seen from Figure 
5-24, the presence of excess number of the numeral 3 confirms the results obtained 
from the curvature maps for the sphere range data. Figure 5-25 is the plot using the 
coefficients listed in Table 5-2 for the cylinder. Both of these plots validate our 
findings from the analysis of the curvature maps. 

The experiments mentioned above were performed on a large number of real 
range data sets for spheres, cylinders, and cones, the results of which are shown in 
appendix A. 

5.3 Application of the Recognition Process to the Processed Image Data 

The next objective was to apply our recognition schemes to the processed images 
of a sphere, a cylinder and a quadric cone. Each of the processed images of the 
sphere, the cylinder, and the cone were intersected with two planes (one parallel to the 
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Figure 5-24. Best fit plot for the sphere raw image. 
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Figure 5-25. Best fit plot obtained for the cylinder belonging to set A. Numerals 
"1, 3, 5, 7" denote the original image, the 3 x 3 image, the 5 x 5 image, and the 
7x7 image respectively. 
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COMPARISON OF COEFFICIENTS EVALUATED FOR THE ORIGINAL 

AND THE PROCESSED IMAGES OF A SPHERE 

COEFFICIENT 

RAW IMAGE 

3X3 FTLT. IMAGE 

5X5 FILT. IMAGE 

7X7 FILT. IMAGE 

A. COEFF. OF X 2 

0.3026 

0.2211 

-0.3260 

0.4242 

B, COEFF. OF Y 2 

0.2734 

0.2802 

-0.4860 

0.2178 

C, COEFF. OF Z 2 

0.6545 

0.7747 

-0.3338 

0.5845 

F, COEFF. OF YZ 

0.5310 

-0.5348 

0.4834 

-0.3417 

G, COEFF. OF XZ 

0.6357 

-0.4860 

0.7194 

-0.7452 

H, COEFF. OF XY 

0.3524 

0.2339 

-0.5801 

0.4353 

P, COEFF. OF X 

0.3036 

0.1999 

-0.3159 

0.3127 

Q, COEFF. OF Y 

0.4199 

0.4401 

-0.3524 

0.1996 

R, COEFF. OF Z 

-0.8172 

-1.0163 

0.3191 

-0.5858 

D, CONSTANT 

0.2847 

0.3717 

-0.0973 

0.1516 


Table 5-1. Comparison of the coefficients evaluated for the original and the processed 
images of a sphere. 
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COMPARISON OF COEFFICIENTS EVALUATED FOR THE ORIGINAL 

AND THE PROCESSED IMAGES OF A CYLINDER 

COEFFICIENT 

RAW IMAGE 

3X3 FILTERED IMAGE 

5X5 FILTERED IMAGE 

A. COEFF. OF X 1 

0.8338 

0.6636 

0.0572 

B. COEFF. OF Y ! 

0.0041 

0.0209 

0.599 

C. COEFF. OF Z ! 

0.059 

-0.0923 

0.4416 

F, COEFF. OF YZ 

-0.00103 

-0.0219 

-0.807 

G, COEFF. OF XZ 

-0.636 

-0.7604 

0.459 

H, COEFF. OF XY 

0.4437 

0.7727 

-0.149 

P. COEFF. OF X 

-0.0141 

0.4242 

-0.5915 

Q, COEFF. OF Y 

-0.189 

-0.2155 

1.089 

R, COEFF. OF Z 

0.193 

0.374 

-1.019 

D, CONSTANT 

-0.1341 

-0.253 

0.664 


Table 5-2. Comparison of the coefficients evaluated for the original and the 
processed images of a cylinder. 
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xy plane, the other parallel to the xz plane). The results obtained for the sphere are 
tabulated in Table 5-3. A decision on the curve being an ellipse or a circle was made 
based upon the parity and disparity of the x 2 , y 2 , and z 2 coefficients. 


Table 5-3. Curves intercepted by the two planes, z = k, y = k, with real raw and 
processed range data of the sphere. 

Experiments conducted with the raw and the processed images of the cylinder led 
to the results tabulated in Table 5-4. 


Table 5-4. Curves intercepted by the two planes, z = k, y - k, with real raw and 
processed range data of the cylinder. 

As seen from the tabulated results, the raw images come close in generating the 
desired curves for each of the objects, but at the same time a 5 x 5 filter in either case 
generates the exact two-dimensional curves. 

Experiments conducted with the raw and the processed images of the quadric 
cone led to the results tabulated in Table 5-5. 


Cylinder Images 

plane 1, y = k I Ellipse I Line I Line 

plane 2, z = k Ellipse Ellipse Ellipse 


Sphere Images 

plane l,y = k I Ellipse I Ellipse I Ellipse 
plane 2, z = k Ellipse Ellipse Circle 
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Quadric cone Images 

plane 1, z = k 

Ellipse 

Ellipse 

Ellipse 

plane 2, y = k 

Ellipse 

Hyperbola 

Hyperbola 


Table 5-5. Curves intercepted by the two planes, z = k, y - k, with real raw and 
processed range data of the quadric cone. 

As seen from the tabulated results, the raw images do not come close in generating the 
desired curves for each of the objects, but at the same time the 3x3 as well as the 
5x5 filter in either case generates the exact two-dimensional curves. Coefficients 
generated for the raw image data of the cone as well as the 3 x 3 and 5x5 median 
filtered image data of the quadric cone are listed in Appendix B. 

5.4 Application of the Rotation Alignment Algorithm 

The rotation alignment algorithm which determines the orientation of the quadric 
surfaces in space and then aligns them in accordance to our desired coordinate system 
was applied to simulated data as well as real data. 

Consider tables 5-6, 5-7, and 5-8, which compare the coefficients of the sphere 
range data before and after rotation alignment. 

Each of the image data sets, i.e., the original raw image of the sphere and the 3 
x 3 and the 5x5 median filtered images of the sphere, required three iterations to 
eliminate the product terms. Since a sphere is symmetric about all coordinate axes, no 
rotation alignment should be needed. The alignment algorithm was performed just to 
see how the coefficients relate to each other before and after the rotation. The 
coefficients were basically invariant as expected. 

Similarly, tables 5-9, 5-10, and 5-11, show the coefficients obtained before and 
after the rotation alignment for the cylinder range data. However, in these cases, each 
of the image data sets, i.e., the original raw image of the cylinder, the 3x3 and the 
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_ 

COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 


A, COEFF. OF X 2 

0.3026 

0.3206466 

— 

B, COEFF. OF Y 2 

0.2734 

0.184263 


C, COEFF. OF Z 2 

0.6545 

0.7999953 


F. COEFF. OF YZ 

0.5310 

0.0 


G. COEFF. OF XZ 

0.6357 

0.0 


H, COEFF. OF XY 

0.3524 

0.0 

— 

P, COEFF. OF X 

0.3036 

0.252 


Q, COEFF. OF Y 

0.4199 

0.41686 

- 

R, COEFF. OF Z 

-0.8172 

-0.8623 


D, CONSTANT 

0.2847 

0.2847 


Table 5-6. New coefficients of the raw image data of sphere after alignment. 


COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A. COEFF. OF X 2 

0.264 

0.249 

B, COEFF. OF Y 2 

0.129 ! 

0.1311 

C, COEFF. OF Z 2 

0.5738 

0.634 

F, COEFF. OF YZ 

-0.6275 

0.0 

G, COEFF. OF XZ 

-0.783 

0.0 

H, COEFF. OF XY 

0.4014 

0.0 

P, COEFF. OF X 

0.4826 

0.4405 

Q, COEFF. OF Y 

0.3670 

0.3746 

R. COEFF. OF Z 

-0.7218 

-0.7401 

D, CONSTANT 

0.2210 

0.2210 


Table 5-7. New coefficients of the 3 x 3 median filtered image data of sphere 
after alignment. 
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COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.303 

0.3765 

B, COEFF. OF Y 2 

0.392 

0.372 

C, COEFF. OF Z 2 

0.6526 

0.6417 

F, COEFF. OF YZ 

-0.4487 

0.0 

G. COEFF. OF XZ 

-0.8376 

0.0 

H, COEFF. OF XY 

0.2416 

0.0 

P, COEFF. OF X 

0.4047 

0.4259 

1 

Q. COEFF. OF Y 

0.2214 

0.2184 

R, COEFF. OF Z 

-0.7089 

-0.7423 

D, CONSTANT 

0.1913 

0.1913 


Table 5-8. New coefficients of the 5 x 5 median filtered image data of sphere 
after alignment. 


COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.8338 

0.9812 

B, COEFF. OF Y 2 

0.00411 

-0.1143 

C, COEFF. OF Z 2 

0.059 

0.03255 

F, COEFF. OF YZ 

-0.00103 

0.0 

G, COEFF. OF XZ 

-0.636 

0.0 

H, COEFF. OF XY 

0.4437 

0.0 

P, COEFF. OF X 

-0.0141 

-0.0876 

Q, COEFF. OF Y 

-0.189 

-0.2478 

R, COEFF. OF Z 

0.193 

-0.01125 

D. CONSTANT 

-0.1341 

-0.1341 


Table 5-9. New coefficients of the raw image data of cylinder after alignment. 
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COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.6636 

0.9235 

B, COEFF. OF Y 2 

0.0209 

- 0.2967 

C, COEFF. OF Z 2 

-0.0923 

-0.00207 

F, COEFF. OF YZ 

-0.0219 

0.0 

G, COEFF. OF XZ 

-0.7604 

0.0 

H, COEFF. OF XY 

0.7727 

0.0 

P, COEFF. OF X 

0.4242 

0.1923 

Q, COEFF. OF Y 

-0.2155 

-0.5368 

R, COEFF. OF Z 

0.374 

0.05379 

D, CONSTANT 

-0.253 

-0.2533 


Table 5-10. New coefficients of the 3 x 3 median filtered image data of cylinder 
after alignment. 


COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.0572 

-0.07251 

B, COEFF. OF Y 2 

0.599 

0.977 

C. COEFF. OF Z 2 

0.4416 

0.1930 

F, COEFF. OF YZ 

-0.807 

0.0 

G, COEFF. OF XZ 

0.459 

0.0 

H, COEFF. OF XY 

-0.149 

0.0 

P, COEFF. OF X 

-0.5915 

-0.1764 

Q, COEFF. OF Y 

1.089 

1.5696 

R. COEFF. OF Z 

-1.019 

-0.1902 

D, CONSTANT 

0.664 

0.664 


Table 5-11. New coefficients of the 5 x 5 median filtered image data of cylinder 
after alignment. 
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5x5 median filtered images of the cylinder, required four iterations to eliminate the 
product terms. Except for the coefficients of the raw image data, both the 3x3 and 
5x5 filtered coefficients after alignment yielded the desired curves when intersected 
with various planes. Making use of the coefficients of the 5x5 filtered image, the 
diameter of this particular cylinder was calculated to be 4.99 centimeters. The actual 
diameter of the cylinder was 5 centimeters. Appendix C contains more results 
obtained while carrying out the rotation alignment algorithm for other cylinder range 
images. 

The rotation alignment technique was utilized for a large group of simulated data. 
Listed in tables 5-12, 5-13, and 5-14 are several upon which the utilization of our 
recognition scheme correctly identified the surfaces. Upon application of our recogni- 
tion scheme the quadric surfaces represented in Tables 5-12, 5-13, and 5-14 were 
correctly recognized as an ellipsoid, a hyperboloid of one sheet, and a hyperbolic 
cylinder, respectively. 

All the simulated data sets of quadric surfaces could be recognized after conduct- 
ing the rotation alignment technique on the original quadratic representation. The 
three-dimensional discriminant approach which was described in Chapter Two was 
applied to several simulated data of quadrics. 

5.5 Application of Three-Dimensional Discriminant Technique 

Results for the simulated data are illustrated in Table 5-15 and are very effective 
as predicted by the theory. Object (1) refers to a parabolic cylinder, (2) refers to a 
hyperbolic paraboloid, (3) refers to a hyperboloid of one sheet, (4) refers to an ellip- 
soid, (5) refers to a hyperbolic cylinder, (6) refers to a quadric cone, (7) refers to a 
hyperboloid of two sheets, and (8) refers to an elliptic paraboloid. A listing of a sam- 
ple data file is included in Appendix D. However, as expected, unsatisfactory results 
were obtained while experimenting with real data. A listing of all the programs 
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COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

103 

49.84 

B, COEFF. OF Y 2 

125 

96.887 

C, COEFF. OF Z 2 

66 

145.3905 

F. COEFF. OF YZ 

-60 

0.0 

G, COEFF. OF XZ 

-12 

0.0 

H, COEFF. OF XY 

-48 

0.0 

P, COEFF. OF X 

0.0 

0.0 

Q, COEFF. OF Y 

0.0 

0.0 

R, COEFF. OF Z 

0.0 

0.0 

D, CONSTANT 

-294 

-294 


Table 5-12. New coefficients of an unknown simulated data obtained after alignment. 


COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.0 

2.0 

B, COEFF. OF Y 2 

2.0 

-4.0 

C, COEFF. OF Z 2 

1.0 

-1.0 

F, COEFF. OF YZ 

-4.0 

0.0 

G, COEFF. OF XZ 

-4.0 

0.0 

H, COEFF. OF XY 

0.0 

0.0 

P, COEFF. OF X 

0.0 

0.0 

Q, COEFF. OF Y 

0.0 

0.0 

R. COEFF. OF Z 

0.0 

0.0 

D, CONSTANT 

-4.0 

-4.0 


Table 5-13. New coefficients of an unknown simulated data obtained after alignment. 
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COEFFICIENT 

BEFORE 

AFTER ALIGNMENT 

A, COEFF. OF X 2 

0.0 

3.0 

B, COEFF. OF Y 2 

0.0 

0.0 

C, COEFF. OF Z 2 

0.0 

-3.0 

F, COEFF. OF YZ 

-1.414 

0.0 

G, COEFF. OF XZ 

0.0 

0.0 

H, COEFF. OF XY 

1.0 

0.0 

P, COEFF. OF X 

0.0 

0.0 

Q, COEFF. OF Y 

0.0 

0.0 

R, COEFF. OF Z 

0.0 

0.0 

i 

D. CONSTANT 

-3.0 

-3.0 


Table 5-14. New coefficients of an unknown simulated data obtained after alignment. 
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SURFACE CHARACTERIZATION USING THREE-DIMENSIONAL DISCRIMINANT APPROACH FOR 

SIMULATED DATA 

1 COEFFICIENTS OF THE SIMULATED OBJECTS 



A, COEFF. OF X 2 

1 

0 

0 

l 

0 

0 

1 

3 

B. COEFF. OF Y 1 

4 

1 

0 

3 

0 

0 

0 

0 

C, COEFF. OF Z 1 

9 

20 

0 

2 

6 

1 

3 

2 

F, COEFF. OF YZ 

-6 

-4.5 

1.5 

-1 

1.5 

3 

1 

0 

G, COEFF. OF XZ 

3 

-2.5 

1 

0 

1 

-2 

-1 

-2 

H, COEFF. OF XY 

-2 

0.5 

0.5 

1 

0.5 

1 

0 

0 

P, COEFF. OF X 

1 

0.5 

-1 

2 

-2 

2 

0 

0 

Q, COEFF. OF Y 

7 

0 

0 

5 

3 

3 

1 

1 

R, COEFF. OF Z 

0 

0 

3 

0 

0 

0 

3 

0 

D, CONSTANT 

10 

0 

0 

8 

0 

12 

9 

19 

OBJECT IS 

a) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 


Table 5-15. Surface characterization using three-dimensional discriminant approach 
for simulated data. 
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utilized in this research is included in Appendix E. The program that generates the ten 
coefficients for quadric surfaces and the program which aligns the quadnc surfaces to a 
desired coordinate system are among those listed. 


CHAPTER SIX 


CONCLUSIONS 


6.1 Overview 

We have presented a new approach based on two-dimensional analytic geometry 
to recognize a series of three-dimensional objects. Among the various three- 
dimensional objects considered are the hyperboloids of one and two sheets, the ellip- 
soids, the spheres, the circular and elliptical quadric cones, the circular and elliptical 
cylinders, the parabolic and hyperbolic cylinders, the elliptic and hyperbolic para- 
boloids, and the parallelepipeds. 

Our proposed method utilizes a two-dimensional discriminant which is a measure 
for distinguishing curves. Instead of evaluating the ten generated coefficients and 
attempting to recognize the surface from its quadric representation, we can identify 
the quadrics using the information resulting from the intersection of the surface with 
different planes. If the surface is one of those listed above, there are five possible 
two-dimensional curves that may result from such intersections: (i) a circle, (ii) an 
ellipse, (iii) a parabola, (iv) a hyperbola, and (v) a line. Thus, a feature or pattern 
vector with five independent components can be formed for characterizing each of the 

surfaces. 

Although all of the quadric surfaces considered have been symmetric, our recog- 
nition system can be extended to other three-dimensional objects. Figures 6-1, 6-2, 
and 6-3 are examples of these surfaces which exist in the real world. To recognize 
complex objects a suitable segmentation technique is required for the isolation of each 
individual surface. 
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Figure 6-1. This delta rocket is composed of cylindrical and conical shapes (Courtesy 
NASA ). 
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6.2 Advantages of the Recognition Scheme 

Some of the advantages of our recognition scheme are listed below: 

(1) Recognition systems using the curvature approach (evaluation of the mean and 
Gaussian curvatures) are very computationally intensive. These approaches never 
really describe the quadric surface in question. Our proposed recognition system is 
computationally efficient. All of the quadnc surfaces are recognized as well as 
described in terms of their dimensions. 

(2) The three-dimensional discriminant approach discussed in Chapter Two works 
only on ideal or simulated data. It is not useful for real range data. Our recognition 
system is shown to work for both simulated and real range data. 

(3) The best-fit plot technique used for analyzing processed range images is a new 
and efficient technique to determine the validity and integrity of the processed range 
images. 

(4) The rotation alignment technique is a new method which systematically and 
effectively eliminates the product terms and aligns the quadric surfaces in our desired 
coordinate system through an iterative technique. 

(5) The curvature analysis technique and the best-fit plot can be used to determine 
performances of various laser range mappers. 

The equations of the planes which determine distinct feature vectors for each of 
the quadric surfaces are very sensitive to the quality of the digitized range data. In 
case the coefficient determining algorithm does not perform as expected, errors might 
be encountered while forming the feature sets. Active sensors like laser range mappers 
have only recently been developed. Much improvement is expected in the quality of 
range images in the near future. This will make the various recognition schemes much 


more reliable and flexible. 
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6.3 Future Goals and Research Directions 

As seen from the best-fit plot in Chapter Five, regions within the range images 
have been marked to indicate which particular filter size (median) fits the image data 
the best at a particular pixel. While arriving at the coefficients of a second degree 
polynomial which describe a quadric surface, the range data considered were either the 
original, the 3 x 3, the 5 x 5, or the 7 x 7 median filtered range images. A filter 
whose mask size varies from region to region could possibly be a more effective filter 
which would not significantly distort the images. 

Though experiments were performed on a large sample of range images belong- 
ing to spheres, cylinders, and cones, the effectiveness and accuracy of the developed 
recognition system can be tested further by using real range images of paraboloids, 
hyperboloids and cylinders (hyperbolic and parabolic). The recognition algorithm, 
however, has been very accurate when applied to simulated data. 

We propose to extend our recognition algorithm to recognize quadric surfaces 
from complex scenes (scenes composed of more than two objects). This can be 
achieved by first utilizing an effective (existing) segmentation process, whereby range 
data of various surfaces will be separated. We could extend our recognition system to 
recognize irregular surfaces which are made up piece-by-piece of regular quadric sur- 
faces. 

Finally, we would like to investigate fully the mapping of the extrinsic and intrin- 
sic representations of quadric surfaces. This process will lead to a three-dimensional 
discriminant analogous to the two-dimensional discriminant, which will distinguish all 
of the quadric surfaces considered for the recognition process in the course of this 
research. The development of such a discriminant will not only reduce the computa- 
tional complexity, but will also eliminate the process of eliminating the product terms 
(rotation parameters) from the representation of the quadric surfaces. This approach 
will be invariant to pose and orientation. 
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APPENDIX A 


Appendix A consists of the ten coefficients generated for the original and pro- 
cessed range images of a sphere and a cylinder whose data is mapped using a different 
type of laser range mapper. Files with extension *.cod refer to the range data con- 
verted into cartesian coordinates, and files with extension *.coe consist of the 
coefficients generated for each of the images. 
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The input file was "spavgmedl.cod 
The output file is "spavgmedl.coe 
The coeff of x-squared is 0.2963710 
The coeff of y-squared is -7.1920902E-03 
The coeff of z-squared is 0.6404306 
The coeff of yz is -0.2438449 

The coeff of zx is -0.9575970 

The coeff of xy is 0.1657399 

The coeff of x is 0.5431624 

The coeff of y is 0.1216914 

The coeff of z is -0.6774822 

The constant d is 0. 1 

Coefficients for an averaged sphere image. 


The input file was "SPAMED31.C0D 
The output file is "SPAMED31.C0E 
The coeff of x-squared is 0.1939911 
The coeff of y-squared is -6.4082608E-02 
The coeff of z-squared is 0.7181194 


The coeff of 

yz 

is -0.1474224 

The coeff of 

zx 

is -0.9272834 

The coeff of 

xy 

is 5.9526406E-02 

The coeff of 

X 

is 0.6028971 

The coeff of 

y 

is 6.9603384E-02 

The coeff of 

z 

is -0.8269780 

The constant 

d 

is 0.2249310 


Coefficients for a 3 x 3 median filtered averaged sphere image. 


The input file was "SPAMED51.C0D 
The output file is "SPAMED51.COE 
The coeff of x-squared is 0.2154791 
The coeff of y-squared is 0.1832946 
The coeff of z-squared is 0.6808279 


The coeff of 

yz 

is 0.8267535 

The coeff of 

zx 

is -0.2360811 

The coeff of 

xy 

is -0.4166962 

The coeff of 

X 

is -0.2436151 

The coeff of 

y 

is -0.4302364 

The coeff of 

z 

is -0.4764909 

The constant 

d 

is 9.0547577E-02 


Coefficients for a 5 x 5 median filtered averaged sphere image. 



The input file was "stavgmedl.cod 
The output file is "stavgmedl.coe 
The coeff of x-squared is -0.1682273 
The coeff of y-squared is 4.3279368E-02 
The coeff of z-squared is 0.7034476 


The coeff of 

yz 

is 

4.8151415E-02 

The coeff of 

zx 

is 

0.9669251 

The coeff of 

xy 

is 

0.1127514 

The coeff of 

X 

is 

-0.7436121 

The coeff of 

y 

is 

-8.4567606E-02 

The coeff of 

z 

is 

-1.537530 

The constant 

d 

is 

0.7877931 


Coefficients for an averaged cylinder image. 
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The input file was "STAMED3 1 .COD 
The output file is "STAMED31.COE 
The coeff of x-squared is 0.2759137 
The coeff of y-squared is 2.7527343E-02 
The coeff of z-squared is 0.7029013 
The coeff of yz is 0.1449835 

The coeff of zx is -0.9098228 

The coeff of xy is -9.6383080E-02 

The coeff of x is 0.5634921 

The coeff of y is -8.973 1783E-02 

The coeff of z is -0.8506840 

The constant d is 0.2536311 

Coefficients for a 3 x 3 median filtered averaged cylinder image. 


input file was "STAMED51.C0D 
output file is "STAMED51.COE 
coeff of x-squared is 0.1115851 
coeff of y-squared is 3.1368352E-02 
coeff of z-squared is 0.8936580 


The 

The 

The 

The 

The 


The coeff of 

yz 

The coeff of 

zx 

The coeff of 

xy 

The coeff of 

X 

The coeff of 

y 

The coeff of 

z 

The constant 

d 


is 0.1347357 
is -0.5961419 
is -4.83962 15E-02 
is 0.4117958 
is -9.9320240E-02 
is -1.295335 
is 0.4731036 


Coefficients for a 5 x 5 median filtered averaged cylinder image. 


APPENDIX B 


This appendix consists of the ten coefficients generated for the original and pro- 
cessed range images of a quadric cone. Files with extension *.cod refer to the range 
data converted into cartesian coordinates, and files with extension *.coe consists of the 
coefficients generated for each of the images. 
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The input file was "coner.cod 
The output file is "coner.coe 
The coeff of x-squared is 0.9966836 
The coeff of y-squared is -4.40009 IE-03 
The coeff of z-squared is -1.723930E-03 
The coeff of yz is 2.29927 5E-02 

The coeff of zx is -0.1116501 
The coeff of xy is -1.4285150E-02 

The coeff of x is 9.4580045E-04 

The coeff of y is -4.749467 6E-03 

The coeff of z is 1.7082826E-03 

The constant d is -2.6372296E-04 

Coefficients for a raw quadric cone image. 
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The input file was "conep.cod 
The output file is "conep.coe 
The coeff of x-squared is 0.9950956 
The coeff of y-squared is -3.4555 167E-02 
The coeff of z-squared is -8.4933 117E-03 
The coeff of yz is 5.0487362E-02 

The coeff of zx is -0.1104977 

The coeff of xy is -4.7736488E-02 

The coeff of x is 9.5897805E-04 

The coeff of y is -1.6880523E-02 

The coeff of z is 6.8076607E-03 

The constant d is -1. 069648 IE-03 

Coefficients for a median filtered quadric cone image. 


APPENDIX C 


This appendix consists of a sample executed file generated using the surface 
alignment algorithm. The coefficients considered are that of a 3 x 3 filtered image of 
a raw cylinder. 
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OUTPUT DATA FILE OF THE SURFACE ALIGNMENT PROGRAM 


****4 


THE COEFFICIENTS CONSIDERED ARE OF THE 3X3 FILTERED IMAGE OF THE RAW 

*^^?*^****************************** M **************************** 

THE NUMBER OF ITERATIONS COMPLETED IS : 3 

*********************** 


COEFF. OF X SQUARE TERM IS : -0.5819000 
COEFF. OF Y SQUARE TERM IS : -2.5060000E-02 
COEFF. OF Z SQUARE TERM IS : -0.4078000 
COEFF. OF YZ SQUARE TERM IS : -9.1289997E-02 
COEFF. OF XZ SQUARE TERM IS : 0.9860000 
COEFF. OF XY SQUARE TERM IS : 8.9539997E-02 
COEFF. OF X TERM IS : -0.3951000 
COEFF. OF Y TERM IS : 4.5000002E-02 
COEFF. OF Z TERM IS : 0.3026000 
CONSTANT OF PROP. IS : -4.48 1000 IE-02 

**************************************************** 

**************************************************** 

**************************************************** 


************ 

************ 

************ 


-0.5854765 -2.1506138E-02 -0.4078000 - 1.248 1645E-02 0.9838626 

00000000E+O0 -0.3974287 1.3208531E-02 0.3026000 

-0 5854765 - 2 . 1405339E-02 -0.4079008 O.OOOOOOOE+OO 0.9837343 


-1.5888708E-02 -0.3974287 
-0.9965052 -2.1405339E-02 

-1.2192142E-02 -0.4991140 
-0.9965433 -2.1367228E-02 

0.0000000E+00 -0.4990523 
-0.9965433 -2.2384373E-02 

-1.2471 168E-05 -0.4990523 
-0.9965433 -2.2384373E-02 

-1.2471168E-05 -0.4990530 
-0.9965433 -2.2384373E-02 

O.OOOOOOOE+OO -0.4990530 
-0.9965433 -2.2384373E-02 

-1.8273728E-23 -0.4990530 


8.3200261E-03 0.3027738 
3.1279027E-03 -1.0188361E-02 0.0000000E+00 
8.3200261E-03 -2.2511929E-02 
3.1279027E-03 -1.0188161E-02 -6.369 1252E-05 
1.1 4400 12E-02 -2.2511929E-02 
4.1450467E-03 0.0000000E+00 -6.2458355E-05 
6.81 0578 IE-03 -2.4316186E-02 
4.1450476E-03 3.8919640E-10 O.OOOOOOOE+OO 
6.8105781E-03 -2.43006 13E-02 
4.1450476E-03 3.8919640E-10 2.49 1243 IE- 15 
6.8137725E-03 -2.4300613E-02 
4.1450476E-03 O.OOOOOOOE+OO 2.4912431E-15 


6.8137725E-03 -2.4300613E-02 
0.0000000E+00 O.OOOOOOOE+OO O.OOOOOOOE+OO O.OOOOOOOE+OO O.OOOOOOOE+OO 
00000000E+O0 O.OOOOOOOE+OO 0.0000000E+00 O.OOOOOOOE+OO 
‘ **************************************************** ************ 
**************************************************** ************ 
**************************************************** ************ 


THE NEW COEFF. OF X SQUARE TERM IS : -0.9965433 
THE NEW COEFF. OF Y SQUARE TERM IS : -2.2384373E-02 
THE NEW COEFF. OF Z SQUARE TERM IS : 4.1450476E-03 
THE NEW COEFF. OF X TERM IS : -0.4990530 
THE NEW COEFF. OF Y TERM IS : 6.8137725E-03 

THE NEW COEFF. OF Z TERM IS : -2.4300613E-02 
THE NEW CONSTANT OF PROP. IS : -4 .48 1000 IE-02 

**************************************************** 


************ 
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**************************************************** ************ 
ABCFGHPQ 


-0.01248 
0.00000 
-0.01019 
-0.01019 
0.00000 
0.00000 
0.00000 
0.00000 

0.00000 

****************************** ********************** 
**************************************************** 


-0.58548 

-0.58548 

-0.99651 

-0.99654 

-0.99654 

-0.99654 

-0.99654 

-0.99654 

0.00000 


-0.02151 

-0.02141 

-0.02141 

-0.02137 

-0.02238 

-0.02238 

-0.02238 

-0.02238 

0.00000 


-0.40780 

-0.40790 

0.00313 

0.00313 

0.00415 

0.00415 

0.00415 

0.00415 

0.00000 


0.98386 

0.98373 

0.00000 

-0.00006 

-0.00006 

0.00000 

0.00000 

0.00000 

0.00000 


0.00000 

-0.01589 

-0.01219 

0.00000 

- 0.00001 

- 0.00001 

0.00000 

0.00000 

0.00000 


-0.39743 

-0.39743 

-0.49911 

-0.49905 

-0.49905 

-0.49905 

-0.49905 

-0.49905 

0.00000 


0.01321 
0.00832 
0.00832 
0.01144 
0.00681 
0.00681 
0.00681 
0.00681 
0.00000 
************ 
************ 


0.30260 

0.30277 

-0.02251 

-0.02251 

-0.02432 

-0.02430 

-0.02430 

-0.02430 

0.00000 


**************************************************** ************ 
**************************************************** ************ 
**************************************************** ************ 


Alpha 


Beta 


Gamma 


4.567488 0.9253277 39.88381 

-0.3581797 -11.29185 -1.7880693E-03 

-3.6674985E-04 4.202751 IE-07 O.OOOOOOOE+OO 

********************************** ****************** 


************ 

************ 


ajlptot bettot gamtot 

4.208942 -10.36652 39.88202 

****** ******************************************** ** 
********************************** ****** ************ 
**************************************************** 
**************************************************** 


************ 

************ 

************ 

************ 


THE ROTATION MATRIX IS : 

0.7640848 -7.1428910E-02 -0.641 1492 

7.9622917E-02 0.9966942 -1. 61493 14E-02 

0.6401832 -3.8710725E-02 0.7672463 



APPENDIX D 


This appendix consists of a sample data file which is generated while executing 
the 3-D discriminant algorithm. The unknown quadric surface is later classified as a 
parabolic cylinder. 
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£***********»** SAMPLE DATA OF 3-D DISCRIMINANT PROGRAM 


Coeff. of x A 2 (A): 
A = 


1 


Coeff. of y A 2 (B): 
B = 


4 


Coeff. of z A 2 (C): 
C = 


9 


Coeff. of yz (F): 
F = 


-6 


Coeff. of xz (G): 
G = 


3 


Coeff. of xy (H): 
H = 

-2 


Coeff. of x (P): 
P = 


1 


Coeff. of y (Q): 
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Q = 


7 


Coeff. of z (R): 
R = 


0 


Constant of prop. (D): 


10 


1 -2 3 

-2 4 -6 

3-6 9 


EE = 

1-231 
-2 4-6 7 

3-690 
1 7 0 10 


dt_e = 
0 


dt_EE = 


0 
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K_K = 

- 0.0000 

0.0000 

14.0000 

rho_3 = 


1 

rho_4 = 

3 

s_d_EE = 
0 

si = 

-1 

s2 = 


1 

s3 = 

1 

flag = 


0 


The sign of the ch. roots are not the same 


The rank of EE is : 1.0000 

The rank of e is : 3.0000 

The sign of the determinant of EE is : 0.0000 

The characterstics roots have the same sign? : 
The object is a PARABOLIC CYLINDER 


0.0000 


APPENDIX E 


This appendix consists of the listings of the following programs: 

1. Program "Median Filtering", which performs the 3 x 3 and 5 x 5 median filtering 
on range images. 

2. Program "Derivatives" that evaluates the first and second derivative with respect 
to x and y axes of the data files and then transforms it into a sign map. 

3. Program "Rangediff' that generates the sign map for each of the range images 
based upon the magnitude of the range value of neighboring pixels. 

4. Program "Surface" that generates the ten coefficients which describe each of the 
range images. 

5. Program "Surface Alignment" which eliminates the product terms from the 
representation of quadric surfaces thereby aligning them according to a desired 
coordinate system. 

6. Program "3-D discriminant" which implements the classification of the quadrics 
based upon the discriminant procedure. 
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C PROGRAM MEDIAN FILTERING 

PARAMETER (N=512) 

INTEGER*2 A(N,N),MED(N,N) 

CHARACTER* 12 INFILE,OUTFILE 
C 

C MAIN PROGRAM 

C 

WRITE(M23) 

123 FORMAT(5X,’INPUT FILE NAME : INFILE’) 
READ(**)INFILE 
WRITE(*,223) 

223 FORMAT(5X,’OUTPUT FILENAME : OUTFILE’) 
READ(*,*)OUTFILE 


OPEN (UNIT = 1 ,FILE=INFILE,RECL=2048 ,ST ATU S = ’ OLD ’ ) 
READ (1,9)((A(I,J),J=1,N),I=1,N) 

9 FORM AT(5 1214) 

M=3 

c CLOSE(l,DISPOSE=’SAVE’) 

CALL MEDFLT(A,MED,N,M) 

OPEN (UNIT =2,FILE=OUTFILE,RECL=2048 ,ST ATUS= ’ NEW ’) 

WRITE (2,11)((MED(I,J),J=1.N),I=1,N) 

11 FORMAT(512I4) 

c CLOSE(2,DISPOSE=’SAVE’) 

STOP 

END 

CC 

CC SUBROUTINE MEDIAN FILTER 

CC 

SUBROUTINE MEDFLT(A,MED,N,M) 

INTEGER*2 A(N,N),MED(N,N),SORT(50) 

LOGICAL NEXCHAN 
C 
C 
C 

MM=M ** 2 

X=(M+l)/2 

Y=X-1 

Ml=(MM+l)/2 
DO 7 I=X,(N-Y) 

DO 9 J=X,(N-Y) 
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K1=0 

DO 11 K=(I-Y),(I+Y) 

DO 13 L=(J-Y),(J+Y) 

K1=K1+1 

SORT(Kl)=A(K,L) 

13 CONTINUE 
11 CONTINUE 

DO 15 I1=1,(MM-1) 

DO 17 K1=1,(MM-I1) 

IF (SORT(Kl).GT.SORT(Kl+l» THEN 
TEMP=SORT(Kl) 
SORT(Kl)=SORT(Kl+l) 
SORT(Kl+l)=TEMP 
END IF 
17 CONTINUE 
15 CONTINUE 

MED(I,J)=SORT(M 1 ) 

9 CONTINUE 

7 CONTINUE 

DO 19 1=1, Y 
DO 21 J=1,N 
MED(U)=A(U) 

MED(N+1-I,J)=A(N+1-I,J) 

MED(J,N+1-I)=A(J,N+1-I) 

MED(J,I)=A(J,I) 

21 CONTINUE 
19 CONTINUE 
RETURN 
END 
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C***** PROGRAM derivatives 

C***** This program determines the derivative along the x and the 
C***** y axes. A group of files can be compared to see whether a 
C***** a particular location has the same curvature or not. 


INTEGER*2 I1,J1,T1,P1,K,L,I,J 

REAL DX 1 ,DX2,DX3,DY 1 ,D Y2,DY3 

REAL DX1 1,DX22,DX33,DY 1 1,DY22,DY33 

REAL D(70,350),E(70,350),A(1000,3).AA(60,50) 

REAL D1(70,350),E1 (70,350) 

INTEGER*2 STREC,ENDREC 

CHARACTER* 12 INFILE 1 ,INFILE2,INFILE3,POINT 
CHARACTER*2 GRAPH1(70,100),GRAPH2(70,100),GRAPH3(70,100) 
CHARACTER*2 GRAPH4(70,100) 

WRITE(*,20) 

20 FORMAT(5X, ’INPUT FILE NAME : INFILE 1’) 

READ(*,*)INFILE1 

OPEN(UNTT=l , FILE=INFILE1, STATUS=’UNKNOWN’) 

DO 100 1=1,969 
READ( 1 ,*)(A(I, J) J =1,3) 

100 CONTINUE 

DO 811 K=l,51 
DO 815 L=l,19 
AA(K,L)=A(L+( 1 9*(K- 1 )),3) 

815 CONTINUE 
811 CONTINUE 

300 FORMAT (5 1 214) 

C** TO FIND THE DERIVATIVE ALONG X-AXIS 
Cllll WRITE(*,908) 

C908 FORMAT(’INPUT THE STARTING RECORD NUMBER: STREC’) 

C READ(*,*)STREC 

C9008 FORMAT(’INPUT THE ENDING RECORD NUMBER: ENDREC’) 

OPEN(UNIT=2,FILE=’FILEl.X’,STATUS=’UNKNOWN’) 

OPEN(UNTT=3,FILE=’FILEl.Y’,STATUS=’UNKNOWN’) 
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OPEN(UNIT=4,FILE=’FILEl.XX’,STATUS=’UNKNOWN’) 

OPEN(UNIT=8,FILE=’FILEl.YY\STATUS=’UNKNOWN’) 

11178 DO 1104 11=1,51 
DO 1204 Jl=l,19 

D(I1J1)=0.5*((AA(I1J1+1)-AA(I1J1))+(AA(I1+1,J1+1)-AA(I1+1J1))) 
D 1 (1 1 , J 1 )=( A A(1 1,J 1 - 1 )-2*(AA(I 1 ,J 1 ))+A A(1 1 ,J 1 +1 )) 


E 1 a U 1 )=( A A(1 1 + 1 ,J 1 )-2* (A A(I 1 ,J 1 ))+A A(1 1- 1 ,J 1 )) 

E(1 1 ,J1)=0.5*((AA(1 1 ,J1 + 1 )- AA(I U 1 +1 ))+(A A(1 1 J 1 )- A A(1 1 + 1 ,J 1 ))) 
1204 CONTINUE 
1 104 CONTINUE 
1965 DO 11104 11=1,51 

WRITE(2, *)(D(1 1 , J 1 ) ,J 1 = 1 , 1 9) 

WRITE(3,*)(E(I 1 ,J 1 ),J 1 = 1 , 1 9) 

WRITE(4, *)(D 1 (1 1 , J 1 ),J 1 = 1 , 1 9) 

WRITE(8,*)(E1 (1 1 ,J 1 ), J 1 = 1 , 1 9) 

11104 CONTINUE 
CLOSE(2) 

CLOSE(3) 

CLOSE(4) 

CLOSE(8) 

OPEN(UNTT=2,FILE=’FILEl.X\STATUS=’UNKNOWN’) 

OPEN(UNIT=3,nLE= , FILEl.Y’,STATUS= , UNKNOWN’) 

OPEN(UNTT=4,FILE=’FILEl.XX\STATUS=’UNKNOWN’) 

OPEN(UNIT=5,nLE=’FILEl.YY’,STATUS=’UNKNOWN’) 

DO 324 11=1,51,1 
READ(2,*)(D(1 1 ,J 1 ),J 1 = 1 , 1 9) 

324 CONTINUE 

DO 325 11=1,51,1 
DO 326 Jl=l,19 

IF (D(I1 ,J 1 ).GT.D(1 1 ,JI+ 1 ))THEN 
GRAPH1(I1,J1)= 

ENDIF 

IF (D(1 1 ,J 1 ).LT.D(1 1 ,JI+ 1 ))THEN 
GRAPH1(I1,J1)= ’+’ 

ENDIF 

IF (D(I1,J1).EQ.D(U,JI+1))THEN 
GRAPH1(U,J1)= ’ ’ 

ENDIF 

326 CONTINUE 

325 CONTINUE 


DO 328 11=1,51,1 
READ(3,*)(D 1(11,11)41=1,19) 

CONTINUE 
DO 329 11=1,51,1 
DO 330 Jl=l,19 

IF (D 1 (1 1 , J 1).GT.D 1 (1 1 , JI+ 1))THEN 
GRAPH2(U Jl)= 

ENDIF 

IF (D 1 (1 1 , J 1 ) . LT. D 1 (1 1 , JI+ 1 ))THEN 
GRAPH2(U,JI)= ’+’ 

ENDIF 

IF (D1(I1,J1).EQ.D1(U,JI+1))THEN 
GRAPH2(U,J1)= ’ ’ 

ENDIF 

CONTINUE 

CONTINUE 

DO 332 11=1,51,1 

READ(4,*)(E(I1,J1),J1=1,19) 

CONTINUE 
DO 333 11=1,51,1 
DO 334 Jl=l ,19 

IF (E(1 1 ,J 1 ) .GT. E(1 1 , JI+ 1 ))THEN 
GRAPH3(U,J1)= 

ENDIF 

IF (E(1 1 ,J 1 ).LT.E(I1 , JI+ 1 ))THEN 

GRAPH3(U,J1)= V 

ENDIF 

IF (E(1 1 ,J 1 ).EQ.E(1 1 ,JI+ 1 ))THEN 
GRAPH3(U Jl)= ’ ’ 

ENDIF 

CONTINUE 

CONTINUE 

DO 336 11=1,51,1 
READ(5 ,*) (E 1 (1 1 ,J 1 ), J 1 = 1 , 1 9) 

CONTINUE 
DO 337 11=1,51,1 
DO 338 Jl=l,19 

IF (E 1 (1 1 ,J 1 ) .GT.E 1 (1 1,JI+ 1 ))THEN 
GRAPH4(I1,J1)= 

ENDIF 

IF (E1(I1,J1).LT.E1(I1,JI+1))THEN 

GRAPH4(U,J1)= ’+’ 

ENDIF 


IF (Elfll,Jl).EQ.El(Il,JI+l))THEN 
GRAPH4(I1,J1)= ’ ’ 

ENDIF 


338 CONTINUE 
337 CONTINUE 
1324 CONTINUE 

OPEN(UNrr=13,FILE=’GRAPH.X’,STATUS=’UNKNOWN ) 
OPEN(UNTT=14,FILE=’GRAPH.Y\STATUS=’UNKNOWN’) 
OPEN(UNn’=15,FILE=’GRAPH.XX , ,STATUS= , UNKNOWN’) 
OPEN(UNrr=16,FILE=’GRAPH.YY’,STATUS=’UNKNOWN’) 

DO 21104 11=1,51,1 

WRITE(1 3 , 1 234)(GRAPH 1 (1 1 ,J 1 ), J 1 = 1 , 19) 
WRITE(14,1234)(GRAPH2(I1,J1),J1=1,19) 

WRITE( 1 5, 1 234)(GRAPH3(1 1 ,J 1 ), J 1 = 1 , 1 9) 

WRITE(1 6, 1 234)(GRAPH4(1 1 ,J 1),J 1 =1 , 19) 

21104 CONTINUE 

1234 FORMAT(30X,20A1) 

C WRITE(*,21) 

C GOTO 64 

END 
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C***** PROGRAM RANGE SIGN MAP (RANGEDIFF) 

£***** jjjjs PROGRAM generates a sign map for data files 

£***** gy TAKING INTO CONSIDERATION THE ABSOLUTE 
C***** DIFFERENCE IN RANGE VALUE OF NEIGHBORING PIXELS. 

INTEGER*2 A(0:511,0:512),D(100,100) 

INTEGERS IlJl.Tl.Pl.ZZ.XX 

CHARACTER* 12 INFILE l,INFILE2,INFILE3,POINT 

CHARACTER*2 GRAPH1(100,100) 

WRITE(*,20) 

20 FORM AT(5X, ’INPUT FILE NAME : INFILE 1’) 

READ(*,*)INFILE1 

OPEN(UNIT=l, FILE=INFILE1, STATUS=’UNKNOWN’, RECL=2048) 

DO 100 1=1,511 

READ( 1 ,300)(A(I,J),J= 1 ,5 12) 

100 CONTINUE 
300 FORMAT(512I4) 

ZZ=1 
C XX=1 

DO 43 1=165,215 
XX=1 

DO 53 J=260,278 
D(ZZ,XX)=A(I,J) 

C ZZ=ZZ+1 

XX=XX+1 
53 CONTINUE 

C XX=1 

ZZ=ZZ+1 
C XX=1 

43 CONTINUE 

WRITE(*,*)XX,ZZ 

OPEN (U NTT =2 ,FI LE= Tan gev al.dat ’ , ST ATU S = ’ UNKN OWN ’ ) 

OPEN(UNIT=3,FILE=Tangediff.dat’,STATUS=’UNKNOWN’) 
c OPEN(UNIT=4,FILE=’FILEl.XX’,STATUS=’UNKNOWN’) 


DO 325 I=1,ZZ-1 
DO 326 J=1,XX-1 
IF (D(I,J).GT.D(I,J+1))THEN 
GRAPH1(I,J)= ’+’ 



ENDIF 

IF (D(I,J).LT.D(I,J+1))THEN 
GRAPH 1(I,J)= 

ENDIF 

IF (D(I,J).EQ.D(I,J+1))THEN 
GRAPH1(IJ)= ’ ’ 

ENDIF 

326 CONTINUE 
325 CONTINUE 

DO 21104 I=1,ZZ-1 

WRITE(3, 1 234)(GRAPH 1 (I,J),J= 1 ,XX- 1 ) 
WRITE(2,3000)(D(I,J),J= 1 ,XX- 1) 

21104 CONTINUE 
1234 FORMAT(35X,20A1) 

3000 FORMAT(I4) 

STOP 

END 
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C Program Surface 

Q ************************************************************** 

C This program approximates the coefficients of a surface 
C generated by given data points. The input file consists of 
C the rectangular coordinates of points on some surface. 
c ************************************************************** 


integer i,j,k,ip 

real x(9000),y(9000),z(9000),x_2(9000) 
real y_2(9000),z_2(9000),p(9000,10) 

real yz(9000),zx(9000),xy(9000),p_ptr(9000,10,10),sc(10,10) 

real a(4,4),b(6,4),b_tr(4,6),c(6,6),h(6,6),h_inv(6,6) 

real ris(4,8),a_inv(4,4),ba_inv(6,4),ba_invbt(6,6),m(6,6) 

real h_invm(6,6),m_pr(6,6),ai(6,6),bi(6,6),ci(6,6) 

real eigval(6,6),eigvec(6,6),ei_vec(6),a_invbt(4,6) 

real alpha(4),beta(6),a_vect(10) 

character* 18 infile.outfile 

Type*,’ Enter coordinates file 

Accept*, infile 

Type*,’ Enter output coefficients file :’ 

Accept*, outfile 

open(unit= 1 ,file=infile,status=’ old ’) 
open(unit=2,file=outfile,status=’new’) 
c ****** The constraint matrix h and hjnv is created ********* 

write(*,3) 

3 format(5x, ’Input total points not exceeding 7750: ip=’) 

read(*,*) ip 
root=l/(sqrt(2.)) 
do 24 i=l,6 
do 26 j=l,6 
h(i,j)=0 
26 continue 
24 continue 
h(l,l)=l 
h(2,2)=l 
h(3,3)=l 
h(4,4)=root 
h(5,5)=root 
h(6,6)=root 
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rootl=sqrt(2.) 
do 20 i=l,6 
do 22 j=l,6 
h_inv(ij)=0 
22 continue 
20 continue 
h_inv(l,l)=l 
h_inv(2,2)=l 
h_inv(3,3)=l 
h_inv(4,4)=rootl 
h_inv(5,5)=rootl 
h_inv(6,6)=rootl 

c****** Data is read here .*.»».*»*******»******•************* 

do 30 i=l,ip 
read(l,*) (x(i),y(i),z(i)) 

30 continue 

c ****** vector P for scatter matrix is formed here ***** 
do 32 i=l,ip 

x_2(i)=x(i)**2 

y_2(i)=y(i)**2 

z_2(i)=z(i)**2 

yz(i)=y(i)*z(i) 

zx(i)=z(i)*x(i) 

xy(i)=x(i)*y(i) 

32 continue 

do 34 i=l,ip 
p(i,l)=x_2(i) 
p(i,2)=y_2(i) 
p(i,3)=z_2(i) 
p(i,4)=yz(i) 
p(i,5)=zx(i) 
p(i,6)=xy(i) 
p(i,7)=x(i) 
p(i,8)=y(i) 
p(i,9)=z(i) 
p(i,10)=l 
34 continue 
do 36 i=l,ip 
do 38 j=l,10 
do 40 k=l,10 

p_ptr(i,j,k)=p(i,j)*p(i,k) 


40 continue 

38 continue 

36 continue 

do 42 j=l,10 
do 44 k=l,10 
sc(j,k)=0 
44 continue 

42 continue 

c **** s ca tter Matrix is formed here ******************* 

do 46 j=l,10 
do 48 k=l,10 
do 50 i=l,ip 

sc(j,k)=sc(j.k)+p_ptr(i,j,k) 

50 continue 

48 continue 

46 continue 

c ******* The Scatter matrix sc is decomposed into a,b,b_tr,c ** 

do 52 i=l,6 
do 54 j=l,6 
c(i,j)=sc(i,j) 

54 continue 

52 continue 

do 56 i=l,6 
do 58 j=l,4 
b(ij)=sc(i,j+6) 

58 continue 

56 continue 

do 60 i=l,4 
do 62 j=l,6 
b_tr(i,j)=sc(i+6,j) 

62 continue 

60 continue 

do 64 i=l,4 
do 66 j=l,4 
a(i,j)=sc(i+6,j+6) 

66 continue 

64 continue 

do 68 i=l,4 
do 70 j=l,4 
ris(ij)=a(i,j) 
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70 continue 

68 continue 

call invers(ris, 4,4,8) 
do 72 i=l,4 
do 74 j=l,4 
a_inv(ij)=ris(i,j) 

74 continue 

72 continue 

c *************** Now to compute M ******************* 
do 76 i=l,6 
do 78 j=l,4 
ba_inv(ij)=0 
78 continue 

76 continue 

do 80 i= 1 ,6 
do 82 j=l,4 
do 84 k=l,4 

ba_inv(i,j)=ba_inv(i,j)+b(i,k)*a_inv(k,j) 

84 continue 

82 continue 

80 continue 

do 86 i=l,6 
do 88 j=l,6 
ba_invbt(i,j)=0 
88 continue 

86 continue 

do 90 i=l,6 
do 92 j=l,6 
do 94 k=l,4 

ba_invbt(i,j)=ba_invbt(i,j)+ba_inv(i,k)*b_tr(k,j) 

94 continue 

92 continue 

90 continue 

do 96 i=l,6 
do 98 j=l,6 

m(i,j)=c(i,j)-ba_invbt(i,j) 

98 continue 

96 continue 

c 

c ******** Now t0 com p U te M’ **************** 
c 

do 100 i=l,6 
do 102 j=l,6 
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h_invm(i,j)=0 
102 continue 

100 continue 

do 104 i=l,6 
do 106 j=l,6 
do 108 k=l,6 

h_invm(i,j)=h_invm(i,j)+h_inv(i,k)*m(k,j) 

108 continue 

106 continue 

104 continue 
do 110 i=l,6 
do 112 j=l,6 


m_pr(i,j)=0 
112 continue 

1 10 continue 

do 114 i=l,6 
do 116 j=l,6 
do 118 k=l,6 

m_pr(i,j)=m_pr(i,j)+h_invm(i,k)*h_inv(k,j) 

118 continue 

116 continue 

114 continue 


c ********* Now to find the eigen values of M’ ********** 
c 

nd=6 

call eig(nd,m_pr,eigval,eigvec) 


c 

q ******* 
^ ******* 
c 


To find the smallest eigen value and its corresponding ** 
eigen vector ***************************************** 


s_eig=eigval(l,l) 
kount=l 
do 120 i=2,6 

if (s_eig.gt.eigval(i,i)) then 
s_eig=eigval(i,i) 
kount=i 
endif 

120 continue 

do 122 1=1.6 

ei_vec(i)=eigvec(i,kount) 
122 continue 
do 124 1=1,6 
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beta(i)=0 
do 126 j=l,6 

beta(i)=beta(i)+h_inv(ij)*ei_vec(j) 

126 continue 
124 continue 
do 128 1=1,4 
do 130 j=l,6 
a_invbt(ij)=0 
do 132 k=l,4 

a_invbt(i,j)=a_invbt(ij)+a_inv(i,k)*b_tr(k,j) 

132 continue 

130 continue 

128 continue 
do 134 1=1,4 
alpha(i)=0 
do 136 j=l,6 

alpha(i)=alpha(i)+a_invbt(i,j)*beta(j) 

136 continue 

alpha(i)=-alpha(i) 

134 continue 

do 138 i=l,6 
a_vect(i)=beta(i) 

138 continue 

do 140 i=l,4 
a_vect(i+6)=alpha(i) 

140 continue 

c do 142 i=l,10 

write(2,*) (’ The input file was "\infile,’"’) 
write(2,*) (’ The output file is "\outfile,””) 
write(2,*) (’ The coeff of x-squared is \a_vect(l)) 
write(2,*) (’ The coeff of y-squared is \a_vect(2)) 
write(2,*) (’ The coeff of z-squared is \a_vect(3)) 
write(2 ,*) (’ The coeff of yz is \a_vect(4» 

write(2,*) (’ The coeff of zx is \a_vect(5)) 

write(2,*) (’ The coeff of xy is \a_vect(6)) 

write(2,*) (’ The coeff of x is \a_vect(7)) 

write(2,*) (’ The coeff of y is ,a_vect(8)) 

write(2,*) (’ The coeff of z is \a_vect(9)) 

write(2,*) (’ The constant d is \a_vect(10)) 
cl42 continue 

close(unit=2,dispose=’save’) 

close(unit=l,dispose=’save’) 

end 


c ********************************************************* 


Subroutine Invers(ris,N,Nx,Mx) 
Dimension ris(Nx,Mx) 

N1=N-1 
N2=2*N 
Do 2 i=l,N 
Do 1 j=l,N 
jl=j+N 

1 ris(ijl)=0. 
jl=i+N 

2 ris(i,jl)=l. 

Do 10 k=l,Nl 

C=ris(k,k) 

If (Abs(C)-O.OOOOO 1 ) 3,3,5 

5 kl=k+l 

Do 6 j=kl,N2 

6 ris(k,j)=ris(k,j)/C 
Do 10 i=kl,N 

C=ris(i,k) 

Do 10 j=kl,N2 
ris(i,j)=ris(i,j)-C*ris(k,j) 

10 Continue 

Npl=N+l 

If (Abs(ris(N,N))-0.00000 1 ) 3,3,19 

19 Do20j=Npl,N2 

20 ris(N,j)=ris(N,j)/ris(N,N) 

Do 200 1=1, N1 

k=N-l 

kl=k+l 

Do 200 i=Npl,N2 
Do 200 j=kl,N 

200 ris(k,i)=ris(k,i)-ris(k,j)*ris(j,i) 

Do 250 i=l,N 
Do 250 j=l,N 
j 1 

250 ris(i,j)=ris(i,j 1) 

Return 

3 Type*,’ Singularity in row found’ 
Return 
End 

Subroutine eig(nd,ai,bi,ci) 
dimension ai(nd,nd),bi(nd,nd),ci(nd,nd) 


integer nl,ml,n2,m2 
nl=nd 
Ml=nd 
n2=nd 
m2=nd 
ANorm=0.0 
Sn=Float(N2) 

Do 100 i=l,N2 
Do 101 j=l,N2 
If (i-j) 72,71,72 

71 Bi(i,j)=1.0 
Goto 101 

72 Bi(i,j)=0.0 
ANorm=ANorm+Ai(i,j)*Ai(i,j) 

101 Continue 

100 Continue 

ANorm=S qrt( A N orm) 

FNorm=ANorm*(1.0E-09/Sn) 

Thr=ANorm 
23 Thr=Thr/Sn 

3 lnd=0 

Do 102 i=2,N2 
il=i-l 

Do 103 j=l,il 

If (Abs(Ai(j,i))-Thr) 103,4,4 

4 Ind=l 
Al=-Ai(j,i) 

Am=(Ai(j,j)-Ai(i,i))/2.0 

Ao=Al/Sqrt((Al*Al)+(Am*Am)) 

If (Am) 5,6,6 

5 Ao=-Ao 

6 Sinx=Ao/Sqrt(2.0*(1.0+Sqrt(1.0-Ao*Ao))) 
Sinx2=Sinx*Sinx 

Cosx=Sqrt(l ,0-Sinx2) 

Cosx2=Cosx*Cosx 
Do 104 k=l,N2 
If (k-j) 7,10,7 

7 If (k-i) 8,10,8 

8 At=Ai(k,j) 

Ai(k,j)=At*Cosx-Ai(k,i)*Sinx 

Ai(k,i)=At*Sinx+Ai(k,i)*Cosx 

10 Bt=Bi(k,j) 

Bi(k,j)=Bt*Cosx-Bi(k,i)*Sinx 

Bi(k,i)=Bt*Sinx+Bi(k,i)*Cosx 
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104 Continue 
Xt=2.0*Ai(j,i)*Sinx*Cosx 
At=Ai(j,j) 

Bt=Ai(i,i) 

Ai(j,j)=At*Cosx2+Bt*Sinx2-Xt 

Ai(i,i)=At*Sinx2+Bt*Cosx2+Xt 

Ai(j,i)=(At-Bt)*Sinx*Cosx+Ai(j,i)*(Cosx2-Sinx2) 

Ai(i,j)=Ai(j,i) 

Do 105 k=l,N2 
Ai(j,k)=Ai(k,j) 

Ai(i,k)=Ai(k,i) 

105 Continue 

103 Continue 

102 Continue 

If (Ind) 20,20,3 

20 If (Thr-FNorm) 25,25,23 
25 Do 110 i=2,N2 

j=i 

29 If ((Abs(Ai(j- 1 ,j- l)))-(Abs(Ai(j,j)») 30,1 10,1 10 

30 At=Ai(j-l,j-l) 

Ai(j-l,j-l)=Ai(j,j) 

Ai(j,j)=At 

Do 111 k=l,N2 


At=Bi(k,j-l) 

Bi(k,j-l)=Bi(kj) 

Bi(k,j)=At 


111 

Continue 


j=j-l 


If (j-1) 110,110,29 

110 

Continue 


do 112 i=l,N2 
do 114 j=l,N2 
ci(ij)=bi(i,j) 
bi(i,j)=ai(i,j) 
1 14 continue 

112 continue 

return 
end 
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c**** PROGRAM SURFACE ALIGNMENT 

( 2 **** This program is used to eliminate the product terms 
C**** from the quadratic representation of any 3D surface. 

C**** The new coefficents generated consisits of the square terms, 
C**** the x, y, z, and the constant term. 


REAL AA,BB,CC,DD,FF,GG,HH,PP,QQ,RR,D,Test_f,Test_g,test_h 
REAL A(50,50),B(50,50),C(50,50),F(50,50) 

REAL G(50,50),H(50,50),ALPHA(100),BETA(100) 

REAL RESULT(200,200),P(50,50),Q(50,50),R(50,50) 

REAL AAA,BBB,CCC,DDD,EEE,FFF,GGG,HHH,III,ROT(3,3) 

REAL DEL1,DEL2,DEL3,A_A,B_B,C_C,F_F,G_G,H_H,GAMMA(100) 
REAL VV,VVV,VWV,VVWV,THRESHLD,INITMIN,ABSA,ABSB,ABSC 

REAL A_AA,B_BB,C_CC,D_DD,P_PP,Q_QQ,R_RR 

REAL ABSF,ABSG,ABSH,ABSP,ABSQ,ABSR,RRR(50),aIptot,bettot 

REAL gamtot 
INTEGER N,M,I,J 

C F(X,Y,Z)=Ax**2+By**2+Cz**2+2Fyz+2Gxz+2Hxy+2Px+2Qy+2Rz+D 

C =0 

C PARAMETER (THRESHLD= 0.00000000000000001) 

OPEN(UNIT= 1 ,FILE=’ CONV ERGENCE.D AT’ ,ST ATU S= ’ NEW ’ ) 

TYPE*, ’ENTER VALUE FOR THRESHLD:’ 

ACCEPT*, THRESHLD 

Type*, ’Enter coef. of x ** 2 :’ 

Accept*,AA 

Type*, ’Enter coef. of y ** 2 
Accept*,BB 

Type*, ’Enter coef. of z ** 2 :’ 

Accept*,CC 

Type*, ’Enter coef. of yz :’ 

Accept* ,FF 

Type*, ’Enter coef. of xz 
Accept*,GG 

Type*, ’Enter coef. of xy :’ 

Accept* ,HH 

Type*, ’Enter coef. of x :’ 

Accept* ,PP 

Type*, ’Enter coef. of y :’ 

Accept* ,QQ 

Type*, ’Enter coef. of z :’ 
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Accept* ,RR 

Type*,’Enter constant of prop. 
Accept* ,D 


A(1,1)=AA 

B(1,1)=BB 

C(1,1)=CC 

F(1,1)=FF 

G(1,1)=GG 

H(1,1)=HH 

P(1,1)=PP 

Q(1,1)=QQ 

R(1,1)=RR 

ABSA=ABS(AA) 

ABSB=ABS(BB) 

ABSC=ABS(CC) 

ABSF=ABS(FF) 

ABSG=ABS(GG) 

ABSH=ABS(HH) 

ABSP=ABS(PP) 

ABSQ=ABS(QQ) 

ABSR=ABS(RR) 

RRR(1)=ABSA 

RRR(2)=ABSB 

RRR(3)=ABSC 

RRR(4)=ABSF 

RRR(5)=ABSG 

RRR(6)=ABSH 

RRR(7)=ABSP 

RRR(8)=ABSQ 

RRR(9)=ABSR 

DO 3980 1=1,9 

IF (RRR(I).EQ.0)THEN 

RRR(I)= 10000 

ENDIF 


3980 CONTINUE 

IN1TMIN=AMIN 1 (RRR( 1 ),RRR(2),RRR(3),RRR(4),RRR(5),RRR(6),RRR(7) 
+ ,RRR(8),RRR(9)) 

WRITE(*,*)INITMIN 
IF (ABS(INTTMIN).LT. 1 .0)THEN 
A(1,1)=A(1,1)/INITMIN 
B(1,1)=B(1,1)/EN1TMIN 


C( 1 , 1 )=C( 1 , 1 )/INITMIN 

F( 1 , 1)=F( 1 , 1 )/INTTMIN 

G(1,1)=G(U)/INITMIN 

H(1,1)=H(1,1)/INITMIN 

P( 1 , 1 )=P( 1 , 1 )/INITMIN 

Q(U)=Q(U)/INITMIN 

Q(1,1)=Q(1,1)/INITMIN 

DD_D=D/INITMIN 

ELSE 

GOTO 3405 
ENDIF 

3405 A(1,1)=AA 

B(1,1)=BB 
C(1,1)=CC 
F(1,1)=FF 
G(1,1)=GG 
H(1,1)=HH 
P(1,1)=PP 
Q(1,1)=QQ 
R(1,1)=RR 


345 if (b(l,l).eq.a(l,l)) then 


c 

c 

c 


57 


goto 1167 
else 
goto 57 
endif 


else 
goto 57 
endif 

alpha(l)=(0.5*ATAND((H(l,l)/(B(l,l)-A(l,l))))) 
A(l,l)=A(l,l)*COSD(ALPHA(l))*COSD(ALPHA(l))+B(l,l) 
SIND(ALPHA(1))*SIND(ALPHA(1))- H(1,1)*SIND(ALPHA(1)) 
COSD(ALPHAU)) 

B(l,l)=B(l,l)*COSD(ALPHA(l))*COSD(ALPHA(l))+A(l,l) 

SIND(ALPHA(1))*SIND(ALPHA(1))+H(1,1)*SIND(ALPHA(1)) 

COSD(ALPHA(l)) 

F(l,l)=G(l'l)*SIND(ALPHA(l))+F(l,l)*COSD(ALPHA(l)) 

G(l,l)=G(U)*COSD(ALPHA(l))-F(l,l)*SIND(ALPHA(l)) 

H(1*1)=0 
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P(1,1)=P(1,1)*C0SD(ALPHA(1))-Q(1,1)*SIND(ALPHA(1)) 

Q(1,1)=Q(1,1)*C0SD(ALPHA(1))+P(1,1)*SIND(ALPHA(1)) 

R(1,1)=R(1,1) 

IF (ABS(F(1,1)).LT.THRESHLD)THEN 

GOTO 1005 

ELSE 

GOTO 1167 
ENDIF 

1005 IF (ABS(G(1,1)).LT.THRESHLD)THEN 
GOTO 1812 
ELSE 

GOTO 1167 
ENDIF 

1167 IF (C(1,1).EQ.B(1,1))THEN 
GOTO 1169 
ELSE 

GOTO 1200 
ENDIF 

1200 BETA(1)=(0.5*ATAND((F(1,1)/(C(1,1)-B(1,1)))» 

A(1,2)=A(1,1) 

B(l,2)=B(l,l)*COSD(BETA(l))*COSD(BETA(l))+C(l,l)* 

+ SIND(BETA(l))*SIND(BETA(l))-F(l,l)*SIND(BETA(l))*COSD(BETA(l)) 


C(l,2)=C(l,l)*COSD(BETA(l))*COSD(BETA(l))+B(l,l)* 

+ SIND(BETA(l))*SIND(BETA(l))+F(l,l)*SIND(BETA(l))*COSD(BETA(l)) 

F(l,2)=0 

G(l,2)=G(l,l)*COSD(BETA(l» 

H(1,2)=-G(1,1)*SIND(BETA(1)) 

P(1,2)=P(1,1) 

Q( 1 ,2)=Q( 1 , 1 )*COSD(BETA( 1 ))-R(l , 1)*SIND(BETA( 1 )) 

R(1 ,2)=R( 1 , 1 )*COSD(BETA(l))+Q(l ,1 )*SIND(BETA( 1 )) 

IF (ABS(H(1,2)).LT.THRESHLD)THEN 

GOTO 1007 

ELSE 

GOTO 1169 
ENDIF 

1007 IF (ABS(G(1,2)).LT.THRESHLD)THEN 
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GOTO 1812 
ELSE 

GOTO 1169 
ENDIF 

1169 IF (C(1,2).EQ.A(1,2))THEN 
GOTO 67 
ELSE 

GOTO 1235 
ENDIF 


1235 GAMMA(1)=(0.5*ATAND((G(1,2)/(C(1,2)-A(1,2))))) 

A(1,3)=A( 1 ,2)*COSD(GAMM A( l))*COSD(G AMMA( 1 ))+C( 1 ,2)* 

+ SIND(GAMMA(1))*SIND(GAMMA(1))-G(1,2)*SIND(GAMMA(1)) 

+ *COSD(GAMMA(l)) 

B(1,3)=B(1,2) 

C(1,3)=C( 1 ,2)*COSD(GAMMA(l))*COSD(GAMMA(l))+A( 1 ,2)* 

+ SIND(GAMMA(1))*SIND(GAMMA(1))+G(1,2)*SIND(GAMMA(1)) 

+ *COSD(GAMMA(l)) 

F(1,3)=H(1,2)*SIND(GAMMA(1)) 

G(l,3)=0 

H(l,3)=H(l,2)*COSD(GAMMA(l)) 

P(l,3)=P(l,2)*COSD(GAMMA(l))-R(l,2)*SIND(GAMMA(l)) 

Q(1,3)=Q(1,2) 

R(1 ,3)=R( 1 ,2)*COSD(GAMMA(l))+P( 1 ,2)*SIND(GAMMA( 1 )) 


IF (ABS(F(1,3)) LT.THRESHLD)THEN 

GOTO 1009 

ELSE 

GOTO 67 

ENDIF 

1009 IF (ABS(H(1,3)).LT.THRESHLD)THEN 
GOTO 1812 
ELSE 
GOTO 67 
ENDIF 
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67 DO 10 1=2,100 

C DO 20 J=1 

71 if((b(i-l,3).eq.a(i-l,3)))then 

goto 167 
C else 

c if(h(i,3).eq.0)then 

c goto 67 

c else 

c goto 67 

c endif 

else 

goto 177 
endif 

177 alpha(I)=(0.5*ATAND((H(I- 1 ,3)/(B(I- 1 ,3)- A(I- 1 ,3))))) 

A(1, 1)=A(I- 1 ,3)*COSD(ALPH A(I))*COSD(ALPH A(I))+(B(I- 1,3))* 
+ SIND(ALPHA(I))*SIND(ALPHA(I))- H(I-1,3)*SIND(ALPHA(I))* 
+ COSD(ALPHA(I)) 

B(I,1)=B(I- 1 ,3)*COSD(ALPHA(I))*COSD(ALPH A(I))+A(I- 1 ,3)* 

+ SIND(ALPHA(I))*SIND(ALPHA(I))+H(I-1,3)*SIND(ALPHA(I))* 
+ COSD(ALPHA(I)) 

ca,i)=ca-i,3) 

F(1, 1 )=F(I- 1 ,3)*COSD(ALPHA(I)) 
G(I,1)=-F(I-1,3)*SIND(ALPHA(I)) 

H(I,1)=0 

P(I,l)=P(I-l,3)*COSD(ALPHA(I))-Q(I-l,3)*SIND(ALPHA(I)) 

Q(I,l)=Q(I-l,3)*COSD(ALPHA(I))+P(I-l,3)*SIND(ALPHA(I)) 

Ra,i)=Ra-i,3) 


IF (ABS(F(I,1)).LT.THRESHLD)THEN 

GOTO 1011 

ELSE 

GOTO 167 

ENDIF 

1011 IF (ABS(G(I,1)).LT.THRESHLD)THEN 
N=I 

GOTO 666 
ELSE 
GOTO 167 
ENDIF 

167 if((c(i,l).eq.b(i,l)))then 
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C 

C 

c 

59 

+ 

+ 


1013 


c69 

69 


goto 69 

else 

goto 59 

endif 

else 

goto 59 

endif 

BETA(I)=(0.5*ATAND((F(I,1)/(C(I,1)-B(I,1))))) 

A(I,2)=A(I,1) 

B(I,2)=B(I,l)*COSD(BETA(I))*COSD(BETA(I))+C(I,l)* 

SIND(BETA(I))*SIND(BETA(I))-F(I,l)*SIND(BETA(I))*COSD(BETA(I)) 

C(I,2)=C(l,l)*COSD(BETA(I))*COSD(BETA(I))+B(I,l)* 

SIND(BETA(I))*SIND(BETA(I))+F(I,l)*SIND(BETA(I))*COSD(BETA(I)) 

F(I,2)=0 

G(I,2)=G(I,l)*COSD(BETA(I)) 

H(I,2)=-G(U)*SIND(BETA(I)) 

P(I,2)=P(I,1) 

Q(I,2)=Q(I,l)*COSD(BETA(I))-Ra.l)*SIND(BETA(I)) 

R(I,2)=R(I,l)*COSD(BETA(l))+Q(I.l)*SIND(BETA(I)) 


IF (ABS(G(I,2)).LT.THRESHLD)THEN 

GOTO 1013 

ELSE 

GOTO 69 

ENDIF 

IF (ABS(H(I,2)).LT.THRESHLD)THEN 
N=I 

GOTO 666 
ELSE 
GOTO 69 
ENDIF 


if(g(i,2).eq.0)then 

if((c(i,2).eq.a(i,2)))then 
goto 10 
else 
goto 63 
endif 


c 

c 
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c else 

goto 63 
endif 

63 GAMMA (I)=(0.5*ATAND((G(I,2)/(C(I,2)-A(I,2))))) 

A(I,3)=A(I,2)*COSD(GAMMA(I))*COSD(GAMMA(I))+C(I,2)* 

+ SIND(GAMMA(I))*SIND(GAMMA(I))-G(I,2)*SIND(GAMMA(I)) 

+ *COSD(GAMMA(I)) 

B(I,3)=B(I,2) 

C(I,3)=C(I,2)*COSD(GAMMA(I))*COSD(GAMMA(I))+A(I,2)* 

+ SIND(GAMMA(I))*SIND(GAMMA(I))+G(U)*SIND(GAMMA(I)) 
+ *COSD(GAMMA(I)) 

F(I,3)=H(I,2)*SIND(GAMMA(I)) 

G(I,3)=0 

H(I,3)=H(I,2)*COSD(GAMMA(I)) 

P(I,3)=P(I,2)*COSD(GAMMA(I))-R(I,2)*SIND(GAMMA(I)) 

Q(I,3)=Q(I,2) 

R(I,3)=R(I,2)*COSD(GAMMA(I))+P(I,2)*SIND(GAMMA(I)) 

IF (ABS(F(I,3)).LT.THRESHLD)THEN 

GOTO 1015 

ELSE 

GOTO 10 

ENDIF 

1015 IF (ABS(H(I,3)).LT.THRESHLD)THEN 
N=I 

GOTO 666 
ELSE 
GOTO 10 
ENDIF 


20 

10 

1812 

666 

123 


CONTINUE 

CONTINUE 

N=1 

WRITE(*,*)N 
WRITER, 123) 

FORMAT(5X ’**************************************************** 
WRrTE(*,*)(’THE NUMBER OF ITERATIONS COMPLETED IS:’,N) 
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M=N*3 

DO 1000 I=1,N 
DO 1001 J=l,3 
RES ULT(3* (I- 1 )+J, 1 )=A(I J) 
RESULT(3*(I-1)+J,2)=B(I,J) 
RESULT(3*(I-1)+J>3)=C(IJ) 
RESULT(3*(I- 1)+J,4)=F(I,J) 
RESULT(3*(I- 1)+J,5)=G(U) 
RESULT(3*(I- 1)+J,6)=H(I,J) 
RESULT(3*(I-1)+J,7)=P(I,J) 
RES ULT(3 * (I- 1 )+ J, 8)=Q(I,J) 
RESULT(3*(I- 1)+J,9)=R(I,J) 


1001 

1000 


CONTINUE 

CONTINUE 

WRITE(1,*)CTHE NUMBER OF ITERATIONS 


WRITE(1,123) 


COMPLETED IS 


N) 


WRITE(1 ,*)(’COEFF. OF X SQUARE TERM IS : AA) 

198 WRITE(l,*)(’COEFF. OF Y SQUARE TERM IS : BB) 

298 WRITE(l,*)(’COEFF. OF Z SQUARE TERM IS : CC) 

398 WRITE(1 ,*)(’COEFF. OF YZ SQUARE TERM IS : \ FF) 

498 WRITE(l,*)CCOEFF. OF XZ SQUARE TERM IS : ’ ,GG) 

598 WRrrE(l,*)(’COEFF. OF XY SQUARE TERM IS : ’ ,HH) 

WRrTE(l,*)( , COEFF. OF X TERM IS : \ PP) 
WRITECl^X’COEFF. OF Y TERM IS : \ QQ) 
WRTTE(l,*)CCOEFF. OF Z TERM IS : RR) 

WRITE(l,*)(’CONSTANT OF PROP. IS : \ D) 

write(l,123) 

write(l,123) 

write(l,123) 

DO 2000 1=1, M 

WRITE(1,*)(RESULT(I,J),J=1,9) 

2000 CONTINUE 

A_AA=RESULT(M-2,1) 

B_BB=RESULT(M-2,2) 

C_CC=RESULT(M-2,3) 

P_PP=RESULT(M-2,7) 
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Q_QQ=RESULT(M-2,8) 

R_RR=RESULT(M-2,9) 

D_DD=D 


do 30001 i=l,3 
write(l,123) 

30001 continue , t ... 

WRITE(1 ,*)(’THE NEW COEFF. OF X SQUARE TERM IS : , A_AA) 

WRITE(1 ,*)(’THE NEW COEFF. OF Y SQUARE TERM IS : \ B_BB) 

WRITE(1,*)(’THE NEW COEFF. OF Z SQUARE TERM IS : ’, C_CC) 

WRITE( 1 ,*)(’THE NEW COEFF. OF X TERM IS : P_PP) 

WRITE(1 ,*)(’THE NEW COEFF. OF Y TERM IS : Q_QQ) 

WRITE(1,*)(’THE NEW COEFF. OF Z TERM IS : R_RR) 

WRITE(1 ,*)(’THE NEW CONSTANT OF PROP. IS : \D_DD) 

do 3001 i— 13 
write(l,123) 

3001 continue 


1278 

+ 

1897 

+ 


write(l,1278) 

format(6x,’A\9x,’B’,9x’C’,9x,’F’,9x,’G ,9x, H ,9x. 


P’, 


9x,’Q\9x,’R’) 
write(l,1897) 
format(5x,’— ■ 


’) 


DO 2001 1=1, M 

WRITE(1,1234)(RESULT(I,J),J=1,9) 

2001 CONTINUE 
1234 format(9F10.5) 

DO 3000 1=1,5 
WRITE(1,123) 

3000 CONTINUE 
write(l,1908) 

1908 format(6x,’Alpha’,9x, ’Beta’, 9x, ’Gamma’) 

write(l,1897) 

DO 4000 1=1, N 

WRITE(1,*)ALPHA(I),BETA(I),GAMMA(I) 
4000 CONTINUE 
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alptot=alpha(l)+alpha(2)+alpha(3) 

bettot=beta ( 1 ) +beta(2) +beta(3) 

gamtot=gamma( 1 )+gamma(2)+gamma(3) 

write(l,123) 

write(l,123) 

write(l,1998) 

1998 format(6x,’ALPTOT\9x,’BETTOT\9x,’GAMTOT) 
write(l ,*)alptot,bettot,gamtot 

write(l,123) 

c***** Yo evaluate coeff. of yz, xz, and xy once alpha, beta 

gamma are evaluated. 

write(*,*)alpha(l),beta(l),gamma(l) 

AAA=BB*cosd(alpha(l))*cosd(alpha(l))+(AA*sind(alpha(l)) 

+ *sind(alpha(l)))+((HH/2)*sind(2*alpha(l)))-CC 
BBB=gg*sind(alpha(l))+(ff*cosd(alpha(l))) 
CCC=((aa-bb))*sind(2*alpha(l))+(hh*cosd(2*alpha(l))) 
DDD=gg*cosd(alpha(l))-(ff*sind(alpha(l))) 

EEE=aa*(cosd(alpha(l ))*cosd(alpha( l))-(sind(alpha(l )) 

+ * sind(alpha( 1 )) * sind(beta( 1 )) * si nd(beta( 1 )))) 

FFF=bb*(sind(alpha( 1 )) * sind(alpha( 1 ))-(cosd(alpha( 1 )) 

+ *cosd(alpha(l))*sind(beta(l))*sind(beta(l)))) 
GGG=cc*cosd(beta(l))*cosd(beta( 1 )) 

HHH=(gg/2)*sind(alpha(l))*sind(2*beta(l))+((ff/2)*cosd(alpha(l)) 

+ sind(2*beta(l))) 

III=(hh/2)*sind(2*alpha(l))*(l+sind(beta(l))* 

+ sind(beta(l))) 

Test_F=(AAA*sind(2*beta(l))+BBB*cosd(2*beta(l)))* 

+ cosd(gamma(l))+(CCC*cos(beta( 1 ))-DDD*sind(beta( 1 ))) 

+ *sind(gamma(l)) 

Test_G=(EEE+FFF-GGG-HHH-III)*SIND(2*GAMMA(l)) + 

+ (CCC*SIND(BETA(l))+COSD(BETA(l))*DDD)*COSD(2*GAMMA(l)) 


TEST H=(CCC*COSD(BETA(l))-DDD*SIND(BETA(l)))*COSD<GAMMA(l)) ; 
+ “(AAA*SIND(2»BETA(l))+BBB*COSD(2*BETA(l)))*SIND(GAMMA(l)) 


c 


write(l,123) 

write(l,123) 
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c write(l,124) 

c 

c write(l ,*)test_f,test_g,test_h 

write(l,123) 
write(l,123) 
write(l,124) 

124 format(5x,’THE ROTATION MATRIX IS : ’) 
c**** To evaluate the rotation matrix 

rot( 1 , 1 )=cosd(alpha( 1 ))*cosd(gamma(l ))-(sind(alpha( 1 ))* 
+ sind(beta( 1 )) * sind(gamma( 1 ))) 

rot(l, 2 )=-sind(alpha(l))*cosd(gamma(l))-(cosd(alpha(l))* 

+ sind(beta(l))*sind(gamma(l))) 

rot( 13 )=-sind(gamrna(l))*cosd(beta(l)) 

rot(2, 1 )=sind(alpha( 1 )) *cosd(beta( 1 )) 
rot( 2 , 2 )=cosd(beta(l))*cosd(alpha(l)) 
rot(2,3)=-sind(beta(l)) 

rot(3, 1 )=cosd(alpha( 1 ))*sind(gamma(l ))+(sind(alpha(l ))* 
+ sind(beta(l))*cosd(gamma(l))) 

rot(3,2)=cosd(alpha(l))*cosd(gamma(l))*sind(beta(l)) 

+ -(sind(alpha( 1 ))*sind(gamma( 1 ))) 

rot(3,3)=cosd(gamma(l))*cosd(beta(l)) 

DO 989 1=1,3 

WRITE(l,*)(ROT(I,J),J=l,3) 

989 CONTINUE 
stop 


end 


C***** PROGRAM 3-D DISCRIMINANT 
C***** Implementation of the 3-D discriminant approach 
C***** implemented on MATLAB 
diary on 

input(’Coeff. of x A 2 (A): ’); 

A=ans 

inputf’Coeff. of y A 2 (B): ’); 

B=ans 

input(’Coeff. of z A 2 (C): ’); 

C=ans 

input(’Coeff. of yz (F): ’); 

F=ans 

input(’Coeff. of xz (G): ’); 

G=ans 

inputfCoeff. of xy (H): ’); 

H=ans 

input(’Coeff. of x (P): ’); 

P=ans 

inputfCoeff. of y (Q): ’); 

Q=ans 

input(’Coeff. of z (R): ’): 

R=ans 

inputfConstant of prop. (D): ’); 

D=ans 

F=F/2; 

G=G/2; 

H=H/2; 

P=P/2; 

0 = 0 / 2 ; 

R=R/2; 

e=[A H G 
H B F 
G F C ] 

EE=[ A H G P 
H B F Q 
G F C R 
P Q R D ] 

dt_e=det(e) 

dt_EE=det(EE) 


K_K=eig(e) 

rho_3=rank(e) 

rho_4=rank(EE) 

s_d_EE=sign(dt_EE) 
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sl=sign(K_K(l)) 
s2=sign(K_K(2)) 
s3=sign(K_K(3)) 
flag=0 
if si = s2; 
flag=flag+l 
end; 

if si = s3; 
flag=flag+l 
end; 

if flag — 2; 
an_w= 1 ; 
fprintf(ViVi 
else; 
an_w=0; 
end; 

fprintf(’ViNn 
fprintf( , \n\n 
fprintf('\n\n 
fprintf(^\n 
fprintf(^Vi 

if rho_3=3 
if rho_4— 4 
if s_d_EE==-l 
if an_w=l 

fprintf(Vi\n The object is an ELLIPSOID \nW) 

end 

end 

end 

end 

if rho_3=3 
if rho_4=4 
if s_d_EE==l 

if an_w=0 ^ v 

fprintff'sn'n The object is a HYPERBOLOID OF ONE SHEET \n\n ) 

end 

end 

end 

end 

if rho_3=3 
if rho_4=4 
if s d_EE==- 1 


if an_w=0 

fprintf(%n\n The object is a HYPERBOLOID OF TWO SHEETS \n ) 

end 

end 

end 


The sign of the ch. roots are the same W) 


The sign of the ch. roots are not the same \n’) 

The rank of EE is : %9.4f \n rho_3 ) 

The rank of e is : %9.4f \n rho_4 ) 

The sign of the determinant of EE is : %9.4f \n\s_d_EE ) 

The characterstics roots have the same sign? : %9.4f \n\ an_w) 
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end 

if rho_3=3 
if rho_4=3 
if an_w=0 

fprintf(’Nn\n The object is a REAL QUADRIC CONE \n ) 

end 

end 

end 

if rho_3=2 
if rho_4=4 
if s_d_EE==-l 
if an_w=0 

fprintf(’\n\n The object is an ELLIPTIC PARABOLOID \n ) 

end 

end 

end 

end 

if rho_3=2 
if rho_4— 4 
if s_d_EE==l 
if an_w=0 

fprintfCNnNn The object is a HYPERBOLIC PARABOLOID Nn ) 

end 

end 

end 

end 

if rho_3=2 
if rho_4=3 
if an_w=l 

fprintf(’\n\n The object is an ELLIPTIC CYLINDER \n ) 

end 

end 

end 

if rho_3=2 
if rho_4=3 

if an_w=0 v 

fprintfC’NnVi The object is a HYPERBOLIC CYLINDER \n ) 

end 

end 

end 


if rho_3=l 
if rho_4=3 

fprintfCNnNn The object is a PARABOLIC CYLINDER \n ) 

end 

end 

diary off 



